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Abstract 


The  configuration  design  space  and  flight  regime  of  new  air  vehicles  is  vast.  Determination  of 
nonlinear  aeroelastic  responses,  vehicle  trim  states  and  vehicle  stability  is  a  challenge  that  must  be 
addressed  in  an  efficient,  methodical  manner.  To  accomplish  this,  we  are  developing  a  method  for 
determining  and  analyzing  nonlinear  aeroelastic  responses  using  a  combination  of  a  continuation 
method  with  a  variable-fidelity  reduced-order  model  of  the  flowfield,  a  nonlinear  description  of 
structural  response,  and  a  trim  model  of  the  air  vehicle  capturing  both  rigid  body  and  flexible 
body  motions. 

This  report  summarizes  the  accomplishments  made  during  this  funding  cycle  of  the  project. 
These  results  include  new  techniques  developed  for  reduced-order  modeling  using  the  proper  orthog¬ 
onal  decomposition  (POD)  method,  a  novel  nonlinear  beam  model,  and  high-quality  grid  generator 
for  aircraft  with  very  large  wing  deformation. 

To  improve  the  efficiency  of  the  reduced-order  model  based  on  POD,  we  developed  several 
acceleration  techniques.  The  most  effective  acceleration  methods  were  database  splitting  and  a  new 
solver  for  linear  systems  with  quasi-symmetric  matrices.  To  improve  the  accuracy  and  efficiency 
of  the  POD  method  applied  to  unsteady  flows  with  shocks,  we  developed  a  method  for  capturing 
discontinuities  using  mathematical  morphology  and  a  technique  for  augmenting  the  POD  basis 
functions. 

To  reduce  the  computational  time  for  structural  modeling  and  to  allow  us  to  use  continuation 
methods,  we  developed  a  novel  nonlinear  beam  model.  Though  some  facets  of  the  method  and  model 
presented  herein  have  been  seen  before,  there  are  several  features  worth  noting  that  distinguish  this 
work.  Compared  to  previous  work,  this  method  offers  improvement  by  accounting  for  a  non-uniform 
beam  with  the  center  of  mass  of  each  cross  section  offset  from  the  elastic  axis.  This  approach 
additionally  permits  rotation  of  the  centroidal  axes  along  the  length  of  the  undeformed  beam 
and  uses  the  Galerkin  method  to  address  abrupt  lengthwise  variations  numerically.  This  model 
improves  upon  prior  work  by  extending  the  order  of  the  nonlinear  terms  retained  to  the  third 
order.  The  method  developed  herein  includes  the  ability  to  account  for  the  nonlinear  contribution 
of  mass-offset,  mass-offset  in  both  dimensions,  non-uniformity,  principal  bending  and  stiffness  axes 
that  vary  mutually  and  along  the  length  of  the  beam,  abrupt  lengthwise  property  changes,  and 
the  consistent  retention  of  all  third-order  terms.  Compared  to  finite  element  results,  the  method 
developed  herein  is  two  orders  of  magnitude  faster,  while  the  frequency  errors  are  typically  less 
than  1%  (the  maximum  error  value  encountered  was  3.4%  for  the  forth  mode). 

Although  the  nonlinear  beam  element  presented  above  excels  on  both  accuracy  and  efficiency, 
not  all  structures  can  be  modeled  by  beams.  For  this  reason,  and  because  we  prefer  to  have  access  to 
the  source  code  rather  than  using  canned  commercial  codes,  during  this  project  we  also  developed 
a  finite  element  method  for  plates.  This  finite  element  code  was  validated  against  NASA  LaRC 
generic  transport  wing  experimental  data. 

To  improve  the  quality  of  the  flow  solver  mesh  generation  for  aircraft  with  extreme  wing  de¬ 
formation  (wing  tip  deformation  up  to  60%  of  the  wing  span)  a  novel  grid  generation  algorithm 
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was  developed  using  conformal  mapping.  Conformal  mapping,  an  almost  extinct  research  topic 
nowadays,  has  the  advantage  of  generating  orthogonal  grids,  that  remain  orthogonal  even  when 
the  wing  deforms.  Consequently,  the  accuracy  of  the  computational  fluid  dynamics  (CFD)  results 
improves.  In  addition,  the  grid  remains  topologically  identical  while  the  wing/fuselage  deform, 
which  simplifies  the  parallel  communication  for  the  flow  solver.  The  results  included  in  this  report 
show  the  improved  grid  quality  of  the  conformal  mapping  grid  compared  to  a  state-of-the-art  grid 
that  uses  Poisson  smoothing. 
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Nomenclature 


Roman  Letters 

A  Cross-sectional  area 
B  Cross  section-fixed  coordinate  system 
B  Constant  in  Weinig-Manea  transformation 
C  Damping  matrix 

C—j  complex  constants  of  conformal  mapping  zi(C) 
c  Cosine 

e v  Offset  of  Centroidal  Axis  above  Elastic  Axis  along  y-axis 
Offset  of  Centroidal  Axis  aft  of  Elastic  Axis  along  ("-axis 
D  Stiffness 
E  Young’s  Modulus 
F  Forcing  Vector 
G  Modulus  of  Rigidity 
j  Cross-Sectional  Mass  Moment  of  Inertia 
K  Stiffness  Matrix 
Torsion  Constant 
L  Length 
£  Lagrangian 
l  Cross-Sectional  Lagrangian 

Number  of  z-direction  Bending  Modes 
M  Mass  Matrix 
m  Cross-Sectional  Mass 

Number  of  y-direction  Bending  Modes 
A f  Inertial  Reference  Coordinate  System 
n  Number  of  Torsional  Modes 
Q  Generalized  Non-Conservative  Sectional  Force 
q  Generalized  Direction 
R  constant  in  Weinig-Manea  transformation 
s  Position  along  Elastic  Axis  of  Beam 
s  Sine 

T  Cross-Sectional  Kinetic  Energy 
T  Kinetic  Energy 
t  Time  or  pitch 
u  Displacement  in  x-direction 

V  Cross-Sectional  Potential  Energy 

V  Potential  Energy 

v  Displacement  in  y-direction 

Vi  y-direction  Displacement  Time  coefficient  for  ith  Out-of-plane  Mode 
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N  V  Inertial  Velocity  of  each  Section 
W  Virtual  Work  Resulting  from  Forces 
Wb  Wnc  at  Boundaries  (s  =  0,  L) 

Wc  Virtual  Work  Resulting  from  Conservative  Forces 
Wi  Mode  Shape  for  the  ith  Bending  Mode 
Wnc  Virtual  Work  Resulting  from  Non- Conservative  Forces 
w  Displacement  in  ^-direction 

Wi  ^-direction  Displacement  Time  coefficient  for  ith  In-Plane  Mode 
x  Inertial  axis  in  lengthwise  direction 
y  Inertial  axis  along  plane  of  fixed  end 

z  Inertial  axis  along  plane  of  fixed  end  or  complex  variable  in  airfoil  domain 
z i  Complex  variable  in  domain  of  the  deformed  circle 
Z2  Complex  variable  in  domain  of  the  deformed  ellipse 


Greek  Letters 

6  Variational  Operator 

(  Cross  section-fixed  axis  or  complex  variable  in  unit  circle  domain 
i]  Cross  Section-Fixed  Axis 
9  Second  Euler  Angle  Rotation  about  yi-axis 
A  Sweep  Angle 
A  Lagrangian  Multiplier 
£  Cross  Section-Fixed  Axis 
p  Mass  Density 
p  Curvature  Vector 

cj)  Third  Euler  Angle  Rotation  about  £-axis 
Torsional  Displacement 

4>i  Torsional  Displacement  Time  Coefficient  for  ith  Torsional  Mode 
ij}  First  Euler  Angle  Rotation  about  z-axis 

a;  Angular  Velocity  of  the  Body  Frame,  Relative  to  the  Inertial  Frame 


Primes  denote  differentiation  with  respect  to  s,  while  over-dots  indicate  differentiation  with  respect 
to  time. 


Chapter  1 

Introduction 


Studies  of  nonlinear  fluid-structure  interactions  have  been  motivated  by  evidence  that  there  are  ad¬ 
verse  aeroelastic  responses  attributed  to  system  nonlinearities.  For  example,  the  flexibility  of  some 
new  classes  of  air  vehicles  leads  to  large  out-of-plane  deformations  as  well  as  to  remarkable  in-plane 
responses  (Strganac  et  al.,  2005).  In  addition,  limit-cycle  oscillations  (LCOs)  occur  in  nonlinear 
aeroelastic  systems  and  remain  a  persistent  problem  on  fighter  aircraft  with  store  configurations. 
Such  LCOs  are  unacceptable  since  aircraft  performance,  aircraft  certification,  mission  capability, 
and  human  factor  issues  such  as  pilot  fatigue  are  adversely  affected.  As  a  result,  high-fidelity  non¬ 
linear  aeroelastic  solvers  for  aircraft  with  large  wing  deformation  are  needed  for  current  and  future 
air  vehicles. 

The  goal  of  this  project  was  to  develop  new  techniques  for  predicting  nonlinear  aeroelastic 
responses.  To  achieve  this  goal  we  developed  novel  methods  for:  (i)  improving  the  accuracy  and 
efficiency  of  reduced-order  modeling  based  on  the  proper  orthogonal  decomposition  (POD)  method, 
(ii)  modeling  nonlinear  structural  dynamics,  and  (iii)  improving  computational  mesh  quality  for 
unsteady  aerodynamics.  This  report  summarizes  these  developments. 

Chapter  2  presents  the  improvements  of  the  reduced-order  model  based  on  the  POD  method. 
The  first  section  describes  a  set  of  acceleration  techniques  that  we  developed  for  the  POD  method. 
The  two  most  efficient  acceleration  techniques  were  included  herein:  database  splitting  and  a  solver 
for  linear  systems  with  quasi-symmetric  matrices.  The  second  section  presents  a  new  POD  method 
for  solving  flows  with  discontinuities,  such  as  shocks. 

Chapter  3  presents  a  novel  nonlinear  beam  model.  The  first  section  describes  the  derivation  of 
the  nonlinear  equations  of  motion,  followed  by  the  modal  representation  which  details  the  deriva¬ 
tion  of  the  linear  and  nonlinear  mass,  stiffness  and  damping  matrices.  Section  three  presents  the 
numerical  implementation  of  the  structural  solver,  followed  by  solver  validation. 

Chapter  4  illustrates  the  use  of  the  continuation  method  for  aeroelasticity  applications.  The 
use  of  the  continuation  method  was  facilitated  by  the  development  of  the  nonlinear  beam  modal 
representation. 

Chapter  5  describes  the  new  grid  generator  for  the  flow  solver.  This  chapter  includes  a  section 
on  the  conformal  mappings  used  to  generate  the  mesh  and  a  section  on  grid  quality.  The  latter 
section  compares  the  grid  quality  of  an  algebraic  mesh  that  was  improved  using  a  Poisson  solver 
with  the  quality  of  the  conformal  mapping  mesh. 
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Chapter  2 

Improvements  of  the  Reduced-Order 
Model  based  on  Proper  Orthogonal 
Decomposition 


This  chapter  presents  the  improvements  done  to  the  POD  method  during  this  project.  The  first 
section  describes  two  acceleration  techniques:  (1)  database  splitting,  and  (2)  a  solver  for  linear 
systems  with  quasi-symmetric  matrices.  The  second  section  presents  a  new  POD  method  for 
solving  flows  with  discontinuities,  such  as  shocks.  This  section  describes  the  method  developed 
for  capturing  discontinuities  using  mathematical  morphology  and  a  new  method  that  we  call  the 
augmented  POD  method.  In  the  augmented  POD  method  the  POD  basis  modes  were  augmented 
by  discontinuity  modes. 


2.1  Acceleration  Techniques  for  the  Proper  Orthogonal  Decom¬ 
position  Method 

Proper  orthogonal  decomposition  (POD)  was  extensively  used  in  aerodynamics  for  extracting  the 
flow  features  from  either  experimental  or  numerically  simulated  data.  The  focus  of  the  work 
presented  herein  was  to  use  the  POD  method  to  develop  a  reduced-order  model  that  diminishes 
the  computational  cost  of  a  full-order  model.  It  should  be  noted  that  when  using  POD,  the 
name  reduced-order  model  is  somewhat  misleading,  because  the  results  produced  by  POD  reduced- 
order  models  are  not  necessarily  lower-fidelity  than  to  those  of  the  full-order  models.  We  have 
demonstrated  that  reduced-order  models  based  on  POD  produce  accurate  results  (shear  stress 
errors  less  than  1%)  for  highly  nonlinear  and  rather  complicated  transport  phenomena,  such  as 
rotor-stator  interaction  in  turbomachinery  (Cizrnas  and  Palacios,  [2003)  and  two-phase  flows  (Yuan 
et  al.,  2005),  if  enough  modes  are  kept  in  the  POD  approximation. 


Existing  POD  implementations  (Lieu  and  Farhat,  2005)  reported  that  the  ROM-based  simula¬ 
tion  was  only  four  times  faster  than  the  full-order  model  simulation.  Lieu  and  Farhat|  (|2005)  also 
reported  that  if  the  simulation  was  limited  to  0.5  seconds,  the  aeroelastic  ROM  was  as  computa¬ 
tionally  intensive  as  the  full-order  nonlinear  aeroelastic  model.  In  this  case,  the  ROM  simulation 
was  beneficial  only  if  longer  time- window  simulations.  We  developed  reduced-order  models  using 
POD  for  viscous  (and  two-phase)  flows  and  obtained  speedups  of  two  orders  of  magnitude  or  more, 
depending  on  the  imposed  accuracy,  that  is,  the  number  of  modes  kept  in  the  solution  ( Cizmas| 
|et  ah,  2008).  It  is  unclear  why  there  is  such  a  large  difference  between  the  speedup  values  reported 
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by  the  Colorado/Stanford  and  Texas  A&M  teams.  It  could  be  that  the  acceleration  techniques 
reported  in  flCizmas  et  al.|,  |200§|)  and  other  practical  implementation  aspects  reported  in  (Premier 
|et  ah,  2009)  are  responsible  for  these  differences. 

The  full-order  model  of  the  unsteady  flow  used  for  the  test  case  utilized  for  developing  accelera¬ 
tion  techniques  consisted  of  the  Reynolds-averaged  Navier-Stokes  equations.  This  full-order  model 
was  used  to  produce  the  POD  basis  functions.  The  energy  spectrum  for  a  channel  flow  with  pres¬ 
sure  variations  of  approximately  10%  of  the  maximum  pressure  value  is  shown  in  Fig.  2T.  It  should 
be  noted  that  the  amount  of  energy  contained  by  the  first  mode  is  rather  large,  varying  between 
48%  to  76%.  These  values  are  larger  than  the  values  obtained  for  other  types  of  flows,  such  as 
turbomachinery  flows  or  two-phase  flows  in  chemical  reactors  where  the  energy  content  of  the  first 
mode  is  less  than  30%.  Consequently,  fewer  modes  are  necessary  to  capture  the  flow  features  for 
aeroelastic  applications,  therefore  leading  to  larger  speed-up  factors. 


Figure  2.1:  Cumulative  energy  spectrum. 

Several  acceleration  methods  of  the  POD  method  have  been  explored.  The  two  most  efficient 
techniques  developed  during  this  project  are  presented  herein:  (i)  an  algorithm  for  splitting  the 
database,  and  (ii)  an  algorithm  for  solving  quasi-symmetrical  matrices.  In  addition  to  these  meth¬ 
ods,  we  also  investigated:  (1)  a  strategy  for  reducing  the  frequency  of  updating  the  matrix  of  the 
system  of  equations  that  generates  the  time  coefficients  of  the  POD,  and  (2)  the  influence  of  the 
time  step  adjustment  strategy  on  the  computational  cost  of  time  integration. 

2.1.1  Database  Splitting 

The  POD  basis  functions  are  extracted  from  a  database  of  snapshots  generated  by  numerically 
integrating  the  governing  differential  equations  of  the  full-order  model.  Currently,  it  is  common 
to  use  a  database  that  includes  all  the  snapshots.  Using  a  single  database  that  covers  the  entire 
time  domain,  however,  could  be  too  restrictive.  For  example,  consider  the  transience  during  the 
start-up  of  the  flow.  The  large  variation  in  time  at  start-up  requires  more  modes  than  are  necessary 
to  model  the  flow  features  present  in  the  latter  part  of  the  simulation.  A  method  to  avoid  this 
problem  is  to  split  the  database  of  snapshots. 
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Splitting  the  database  into  multiple  subsets  produces  an  auto-correlation  matrix  that  contains 
more  relative  energy  in  the  first  modes.  Herein,  energy  (shown  in  Fig.  2.1)  is  defined  as  the 
sum  of  all  the  POD  eigenvalues,  A*.  The  relative  energy  captured  by  the  fcth  mode  is  defined  as 
the  relative  energy  of  the  first  modes  increases,  fewer  POD  modes  are  needed 
in  the  reconstruction  to  approximate  the  solution.  Consequently,  the  computational  cost  of  the 
reduced-order  model  decreases. 

Computing  the  auto-correlation  matrix  for  each  database  subset  is  a  straightforward  process. 
Determining  the  bounds  of  each  subset  to  reduce  the  computational  cost  is  less  trivial.  Two  methods 
for  the  separation  of  the  snapshots  into  subsets  were  explored. 

The  first  method  measured  the  time  variation  of  the  time  coefficients,  a,  of  the  dominant  modes 
of  each  field  variable.  The  variation  of  the  time  coefficients  must  be  calculated  and  monitored  for 
several  modes  for  every  field  variable.  Monitoring  all  of  these  values  can,  in  some  cases,  produce 
conflicting  information.  An  alternative  to  monitoring  all  five  field  variables  was  to  monitor  only 
the  field  variables  that  most  affect  the  flow.  For  the  aeroelastic  case,  the  flow  was  most  affected  by 
the  first  modes  of  velocity  in  the  chordwise  direction. 

The  advantage  of  this  method  is  that  the  end  of  the  transient  regime  can  be  accurately  detected 
during  calculation.  The  disadvantage  of  this  method  is  that  there  is  not  a  unique  value  that 
determines  the  limit  of  the  transience  for  all  five  field  variables. 

The  second  method  proposed  for  separating  the  snapshots  was  to  monitor  the  ratio  of  the  change 
in  CPU  time  and  the  change  in  physical  time,  Atcpu / At.  During  transience,  longer  computation 
times  are  needed  per  time  step  as  shown  in  Fig.  |2.2|.  The  advantage  of  monitoring  this  parameter 
is  that  the  time  slope  is  a  single  value  that  describes  the  behavior  of  all  five  field  variables.  The 
disadvantage  of  this  method  is  that  it  over  predicts  the  end  of  the  transience.  The  magnitude 
of  the  slope  decreased  rapidly  until  t  =  0.3  s  and  continued  to  decrease  somewhat  slower  to  a 
quasi-constant  value  at  t  =  0.45  s.  Placing  the  end  of  the  transient  region  at  t  =  0.45  s  is  more 
conservative  than  the  0.35  s  predicted  by  the  first  method.  Numerical  testing  showed  that  database 
splitting  produced  a  speed-up  factor  of  approximately  1.5. 


Figure  2.2: 


Slope  of  the  total  CPU  time  vs  physical  time  measured  at  0.01  s  increments. 
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2.1.2  Novel  Solver  for  Linear  Systems  with  Quasi-symmetric  Matrices 


Let  us  consider  the  system  of  equations  generated  by  projecting  the  mass  conservation  equation 
onto  the  basis  functions  extracted  from  an  ensemble  of  p  snapshots.  This  system  can  be  written 
as  (pizmas ,  2007) 


Apap  =  & , 

where  ap  is  the  vector  of  unknowns  ap.  The  ik  element  of  the  Ap  matrix  is 

NB 


APlk  =  Wi}T  [A){yk}  -  iAnb]Wknb}, 

nb=  1 


,k  =  l. 


,  m 


(2.1) 


(2.2) 


not  symmetrical  because  of  the  second  term  —  Ylnb=i  {^}T 


[Anbli^knb} 
The  difference  between  the 


The  Ap  matrix  is 

the  element  Apk.  The  matrix  would  be  symmetrical  if  {^fc}  =  { ‘Pk nb] 
two  vectors  {ipk\  and  {<Pknb)  is  small,  however,  because  the  latter  vector  is  evaluated  at  slightly 
different  spatial  locations  compared  to  the  first  vector.  Similarly,  the  systems  of  linear  algebraic 
equations  that  resulted  from  the  momentum  and  energy  conservation  have  matrices  that  are  not 
symmetrical.  These  matrices,  however,  are  quite  close  to  being  symmetrical,  and  for  this  reason 
will  be  called  quasi-symmetrical.  A  typical  example  of  a  quasi-symmetrical  matrix  is  the  A  matrix 
obtained  for  POD  approximation  with  8  modes  (Cizmas,  2007|) 


196.4486 

63.3060 

6.0469 

0.5038 

63.3060 

903.4807 

-44.1690 

6.3410 

6.0459 

-44.1687 

243.2099 

-20.7951 

0.5039 

6.3411 

-20.7953 

930.9194 

-21.3042 

14.0288 

-164.8535 

31.0347 

11.9068 

-7.4940 

68.0527 

20.0167 

2.3477 

6.1634 

19.3267 

14.3861 

-6.8042 

19.8722 

-42.8362 

15.2768 

-21.3047 

11.9071 

2.3488 

-6.8064 

14.0286 

-7.4939 

6.1636 

19.8724 

-164.8536 

68.0529 

19.3275 

-42.8377 

31.0348 

20.0166 

14.3861 

15.2768 

890.8742 

32.1664 

42.8224 

-23.8698 

32.1663 

904.3555 

-10.8230 

26.7999 

42.8222 

-10.8228 

872.6460 

92.5161 

-23.8695 

26.7996 

92.5161 

763.9839 

(2.3) 


The  algorithm  proposed  herein  for  solving  a  system  of  equations  Ax  =  6,  in  which  the  matrix  A 
must  be  positive  definite  and  quasi-symmetrical,  begins  by  splitting  the  matrix  into  a  symmetrical 
and  a  non-symmetrical  part  (pizrnas,  2007): 


{As  +  An)x  =  b 


(2.4) 


The  decomposition  of  the  A  matrix  into  a  symmetrical  and  non-symmetrical  part  is  not  unique. 
This  issue  will  be  discussed  later  in  the  section,  while  for  the  moment  one  will  consider  that  one  of 
the  possible  choices  was  used  to  split  the  A  matrix. 

A  solution  of  the  symmetrical  part  is  then  obtained  by  solving  the  system  of  equations 

A,xW  =  b  (2.5) 

The  solution  of  (p.4|)  is  decomposed  in  a  component  obtained  by  solving  the  system  (|2.5|)  and 
a  correction  needed  because  the  matrix  A  is  non-symmetrical 

x  =  +  x^  (2.6) 
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Substituting  (|2.6D  into  (2A)  and  deducting 

(A,  +  An)x $  =  -An  jW 
(A  +  An)x $  =  b(l) 


yields 


or 


(2.7) 


Note  that  the  system  of  equations  (2.7)  has  the  same  matrix  as  (2.4).  Consequently,  an  iterative 
process  can  be  used  to  find  the  solution,  x^  can  be  considered  a  correction  of  the  solution  xi1'1 
that  is  needed  because  the  A  matrix  is  non-symmetrical.  The  correction  x^  can  be  split  into 
two  components:  a  component  Xg  obtained  by  solving  the  system  AsXg  =  b ^  and  a  correction 
needed  because  the  matrix  A  is  non-symmetrical,  similarly  to  the  approach  used  in  (2.6): 


rU)  =  t(2)  i  t(2) 

Ti  **-'  Q  I  •*-'  rf) 


Substituting  (fO|)  into  (f2.7|)  and  using  Asx^  =  yields 


(A  +  An)xff  —  —A 

(A  +  An)x$  =  fe(2) 


,(2) 


n  ^  s 


or 


This  process  of  approximating  the  solution  yields  after  p  steps 


rn  —  7>(1)  I  ^>(2) 
fXj  |  *Xj  n 


(i)  . 

where  the  values  xy  ,  1  <  i  <  p,  are  obtained  by  solving  the  linear  system 
Asx W  = 


.  T(P)  T(P) 

i  ^  s  i  ? 


(2.8) 


(2.9) 

(2.10) 


(2.11) 


(2.12) 


where  i>W  =  —Anx^  The  iterative  process  of  adding  corrections  is  stopped  when  is  smaller 
than  an  imposed  error. 

The  computation  of  xi^  requires  the  solution  of  the  linear  system  (2.12)  several  times.  The 
Cholesky  decomposition  is  used  for  the  factorization  of  the  positive-definite  symmetrical  matrix 
As.  The  method  takes  advantage  of  the  fact  that  the  As  is  constant  and  only  the  right-hand-side 
vector  iP'1  changes. 

The  method  proposed  herein  replaced  the  LU  decomposition  (press  et  alj,  1992,  p.  34)  that 
was  previously  used  to  solve  the  linear  algebraic  systems  (Yuan  et  al.,  |2005|).  For  a  system  with 
m  equations,  the  number  of  operations  for  the  LU  decomposition  is  m3/ 3  while  the  number  of 
operations  for  the  Cholesky  decomposition  is  m3/ 6. 

As  mentioned  previously  in  this  section,  the  decomposition  of  the  A  matrix  into  a  symmetrical 
and  a  non-symmetrical  part  is  not  unique.  Two  decompositions  were  explored  herein  in  order  to 
evaluate  their  effect  on  the  convergence  of  the  solver  algorithm. 

In  the  first  decomposition,  the  A  matrix,  a,ij,  i,j  =  1, ...  ,m,  was  split  such  that  the  symmetrical 
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and  non-symmetrical  parts  were 


A,  = 


an  a2i 

021  022 

Oml 

Om2 

Oml  Om2  ■  ■ 

Omm 

0  ai2  —  021 

0  0 

•  •  •  Oim 

• • •  O2 m 

Oml 

—  Om2 

0  0 

0 

(2.13) 


An  — 


In  the  second  decomposition,  the  A  matrix  was  split  into  a  symmetrical  matrix  As  and  skew- 
symmetrical  matrix  An 


A*  =  l(A  +  AT) 

An  =  ±(A-AT) 


(2.14) 


These  two  decompositions  were  applied  to  solve  the  system  of  equations  Ax  =  b  where  the  A 
matrix  was  given  by  Q2.3|)  and  the  right-hand-side  term  b  was 


b  =  {-33.601  113.487  -  8.818  39.586  -  21.221  33.225  -  49.507  17.622}T. 

For  an  imposed  error  of  10~9,  the  iterative  solver  algorithm  converged  for  both  decompositions  in 


(2)  .  _ 

three  steps.  The  vector  of  corrections  xy,  shown  in  Table  2.1,  indicates  that  in  these  cases  the 


matrix  decomposition  had  a  minimal  effect  on  the  convergence.  The  Euclidean  norm  of  the  relative 
error  between  the  solutions  obtained  using  the  two  decompositions  was  8.9  •  10-31. 

Numerical  tests  showed  that  three  to  four  terms  in  (2.11)  were  usually  sufficient  to  obtain 


— fi  (i) 

a  solution  with  an  error  less  than  10  b.  Consequently,  three  to  four  xy  solutions  of  the  linear 


algebraic  system  of  equations  (2.12)  must  be  computed.  Since  only  the  right-hand-side  term  lP'> 
changes  while  the  matrix  As  is  constant,  the  number  of  operations  for  the  proposed  method  was 


increased  by  a  factor  proportional  to  m2  multiplied  by  the  number  of  xi*'1  terms  in  (2.11).  As  long 


as  m  is  larger  than  4,  the  computational  cost  of  the  proposed  method  is  approximately  half  that 
of  the  LU  decomposition. 

The  magnitude  of  the  non-symmetrical  terms  was  gradually  increased  in  numerical  tests.  Herein, 


the  degree  of  non-symmetry  was  defined  as  (Li  et  al 
Tj  =  \  \A-  At\\f/ \\ A  +  A1  \\p 


2003|) 


where  ||  \\f  indicates  the  Frobenius  (or  the  Hilbert-Schmidt)  norm  flGolub  and  van  Loan,  1989|, 
p.  56).  The  solution  of  the  Ax  =  b  system  mentioned  above  was  computed  for  three  error  levels: 
10~6, 10— 9 ,  and  10-12.  The  A  matrix  was  decomposed  using  split  option  1  (2.13).  The  degree 
of  non-symmetry  of  the  A  matrix  was  increased  by  multiplying  the  non-symmetric  matrix  As  by 
different  factors  X,  ranging  from  1  to  5- 108.  The  variation  of  the  degree  of  non-symmetry  of  matrix 


A  as  a  function  of  the  factor  X  is  shown  in  Fig.  2.3. 


The  variation  of  the  number  of  iterations  needed  for  convergence  as  a  function  of  the  degree 
of  non-symmetry  77  is  shown  in  Fig.  2.4(a).  It  is  remarkable  that  five  or  less  iterations  are  needed 
for  the  convergence  of  the  solution  for  degree  of  non-symmetry  r/  =  0.01  that  corresponds  to  a 


15 


Table  2.1:  Vector  of  corrections  Xg'1  corresponding  to  matrix  (2.3)  for  decompositions  (|2.13)  and 

(Ill- 


Split  1  (2.13) 

Split  2  (2.14) 

T(i) 

t(2) 

t(3) 

T(i) 

r(2) 

T0) 

s 

1 

0.2205E+00 

-0.5529E-06 

0.7659E-12 

0.2205E+00 

-0.3159E-06 

-0.3772E-11 

2 

-0.1401E+00 

0.3505E-07 

-0.1631E-12 

-0.1401E+00 

0.4948E-07 

0.2289E-12 

3 

0.3053E-01 

-0.2586E-06 

-0.3920E-12 

0.3053E-01 

0.3431E-06 

-0.2669E-11 

4 

-0.4188E-01 

-0.1645E-07 

0.8489E-14 

-0.4188E-01 

0.1406E-07 

0.5404E-14 

5 

0.3669E-01 

-0.8025E-07 

-0.6626E-13 

0.3669E-01 

-0.2367E-07 

-0.6342E-12 

6 

-0.4223E-01 

0.4975E-07 

0.3221E-13 

-0.4223E-01 

0.4024E-07 

0.3478E-12 

7 

0.5685E-01 

0.1344E-07 

0.1147E-13 

0.5685E-01 

0.1726E-06 

0.6542E-13 

8 

-0.1916E-01 

-0.2588E-07 

-0.1567E-13 

-0.1916E-01 

-0.4003E-06 

-0.9619E-13 

multiplying  factor  X  =  10, 000.  The  number  of  iterations  increased  to  approximately  20  for  r/ 
larger  than  0.1  which  correspond  to  multiplying  factors  larger  than  105.  Certainly  at  these  large 
X  values  the  matrix  is  far  from  a  quasi-symmetric  matrix.  If  the  multiplying  factor  X  is  further 
increased,  the  matrix  is  no  longer  positive  definite  and  the  algorithm  cannot  be  used  because  it 
relies  on  the  Cholesky  decomposition. 

The  Eulerian  norm  of  the  difference  between  the  solution  obtained  using  the  LU  decomposition 
and  the  method  proposed  herein  is  shown  in  Fig.  [2.4(b)  as  a  function  of  the  degree  of  non-symmetry 
r/.  The  norm  of  the  difference  increases  as  the  degree  of  non-symmetry  increases.  The  norm  of 
difference  is,  however,  smaller  than  10~10  even  for  values  of  77  as  large  as  0.1.  Consequently,  the 
method  can  be  applied  to  a  larger  class  of  matrices  than  just  matrices  obtained  in  the  proper 
orthogonal  decomposition  method,  as  long  as  these  matrices  are  positive  definite. 


2.2  Proper  Orthogonal  Decomposition  Method  for  Flows  with 
Discontinuities 

2.2.1  Augmented  Proper  Orthogonal  Decomposition  Method 

The  large  computational  time  of  full-order  model  solutions  of  high-fidelity,  unsteady  aerodynamics 
is  a  limiting  factor  for  parametric  studies  in  aeroelasticity.  Note  that  in  aeroelasticity,  the  term 
high-fidelity  aerodynamic  models  currently  means  Reynolds- averaged  Navier-Stokes  solvers.  To 
decrease  computational  time,  reduced-order  models  are  being  developed.  As  illustrated  above,  a 


16 


Figure  2.3:  Degree  of  non-symmetry  rj  as  a  function  of  the  multiplying  factor  X. 


promising  reduced-order  model  is  based  on  proper  orthogonal  decomposition  (POD).  This  section 
presents  a  new  POD  approach  that  can  better  model  flows  with  discontinuities,  such  as  shock 
waves. 

POD  is  a  method  for  extracting  spatially  dependent  orthogonal  basis  functions,  {^(x)}  and 
time  dependent  orthonormal  amplitude  coefficients,  {a(t)},  from  a  function  u(x,  t)  that  varies  in 
both  time  and  space.  The  function  is  usually  parametrized  by  time  to  form  a  collection  of  M 
snapshots,  u(x,tk),k  £  [0,  M].  In  the  POD  approximation, 

m 

u(x,t)  =  ^2aj(tk)<Pj{x),  m  <  M  —  1, 

3=0 


the  basis  functions  are  computed  such  that  the  average  least-squares  truncation  error 

2 


£-m.  — 


«(x,4)  -  ^2aj(tk)<pj(x) 
3=0 


is  a  minimum.  Here  ||  •  |[  denotes  the  L2-norm  given  by  ||/||  =  (,)  denotes  the  Euclidean 

inner  product,  and  (  •  )  denotes  an  ensemble  average  over  the  number  of  observations,  (/)  = 
/(xi  tk)/M  (|Holmes  et  al.,  |l996). 

The  current  POD  methods  produce  unphysical  oscillations  for  flows  with  moving  discontinues. 
These  oscillations  are  due  to  dispersion  errors  caused  by  the  superposition  of  modes  of  increasing 
frequency.  This  phenomenon  is  well  known  and  is  often  referred  to  as  the  Gibbs  phenomenon. 

Lucia  et  al.  ( Lucia  et  al. ,  2003|)  have  studied  the  problem  of  modeling  flows  with  moving  shocks 
using  the  POD  method.  They  have  proposed  a  domain  decomposition  method  where  the  spatial 
domain  was  divided  into  areas  where  shocks  were  possible  and  where  shocks  were  not  possible.  A 
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(a) 


(b) 


Figure  2.4:  Effect  of  degree  of  non-symmetry  on:  (a)  number  of  iterations,  and  (b)  Eulerian  norm 
of  difference  between  the  solutions  of  the  LU  decomposition  and  the  present  method. 


full-order  model  was  solved  where  shocks  were  possible  and  the  rest  of  the  domain  was  solved  using 
a  POD-based  reduced-order  model  (ROM).  This  method  is  somewhat  inefficient,  particularly  when 
a  shock  travels  through  a  large  portion  of  the  spatial  domain. 

To  better  model  flows  with  discontinuities,  we  propose  augmenting  the  POD  database  with 
a  set  of  discontinuity  modes.  The  purpose  of  the  discontinuity  modes  is  to  capture  the  moving 
discontinuities  exactly,  without  the  need  for  domain  decomposition.  The  POD  approximation  can 
then  be  expressed  as 

m  rriA 

u(x,tk)  =  y^Q!j(ti<r)yj(x)  +  y ^Pe(tk)ipe(x-,tk),  (2-15) 

3=0  1=1 


where  'tpe  are  the  discontinuity  modes  and  their  time  coefficients. 

Note  that  the  since  the  POD  method  is  guaranteed  to  produce  orthogonal  basis  functions, 
the  POD  basis  functions  are  linearly  independent.  The  discontinuity  modes,  however,  are  not 
guaranteed  to  be  linearly  independent  of  the  POD  modes.  Therefore,  one  must  verify  the  linear 
independence  of  the  discontinuity  modes  by  appending  them  to  the  POD  modal  matrix  and  com¬ 
puting  the  singular  values.  As  long  as  the  smallest  singular  value  is  larger  than  a  small  number, 
the  appended  modes  are  linearly  independent. 

To  illustrate  the  augmented  POD  approach  we  will  use  the  first-order  wave  equation 


du  du 

~di  +cdi  =  ’ 


(2.16) 


which  is  the  simplest  model  for  a  moving  discontinuity.  This  simple  model  was  chosen  to  test  the 
augmented  POD  method  due  to  the  ease  of  predicting  the  discontinuity  locations 


xs/(U)  =  xs/(t  =  0)  +  cti . 
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Here  the  discontinuity  modes  are  defined  as 


1  x  <  xSte(ti) 
0  x  >  xSte(ti) 


where  xs^{ti)  is  the  location  of  the  7-th  discontinuity  at  time  t  =  ti 
Substituting  (2.15|)  into  (2.16)  yields 


y^dj(4)y>j(x)  +  y ^Pe(tk)ipi(yL,tk)  +  y^^(4)#(x,4) 


j=i 


£=i 


f=i 


+  cya3(4)y'(x)  +  ^PeitkWefatk)  =  0  (2.17) 
j=o 


£=1 


where  ^  :=  §y  and  ft  :=  ^  are  used  for  the  spatial  and  temporal  derivatives,  respectively. 

A  ROM  was  constructed  by  doing  a  Galerkin  projection  of  (|2.17|)  onto  the  augmented  basis 
functions,  which  produced  a  system  of  ordinary  differential  equations  (ODEs) 


lA]U  +[B|o  +M={0}’ 


(2.18) 


where  [A]  <E  Mm+m“  x  {«}  <e  Rm,  {/?}  G  Mm“  and  {d}  G  The  ODE  system 

(|2.18|)  was  solved  using  a  Runge-Kutta-Fehlberg  routine  ( Fchlbcr^,  1969|).  Further  details  of  the 
implementation  are  given  in  (Brenner  et  ah,  2009|) . 

For  this  test  case,  the  initial  velocity  profile  was  a  simple  quadratic  function 

u(x,  t)  =  a  +  b(x  —  ct )  +  d(x  —  ct)2. 

A  single  discontinuity  was  introduced  in  the  initial  condition  according  to 

f  a  +  bx  +  dx2  +  1  x<xsi(t  =  0), 
u[x,  7  =  0)  =  <  ’ 

I  a  +  bx  +  dxz  x  >  xSji(t  =  0). 

Since  only  one  discontinuity  is  present  in  the  flow,  only  one  discontinuity  mode  was  needed,  i.e., 
ma  =  1. 


Figure  2.5  shows  the  velocity  profile  after  one  time  unit  for  a  full-order  model,  a  POD-based 
ROM  with  no  augmentation  and  a  POD-based  ROM  with  augmentation.  The  full-order  model 
solution  was  computed  using  an  explicit  finite  difference  method  that  was  first-order  accurate  in 
space  and  time  and  shows  some  numerical  dissipation  around  the  discontinuity.  The  ROM  without 
augmentation  shows  the  characteristic  oscillations  of  the  Gibbs  phenomenon.  The  ROM  with 
augmentation  produces  the  best  results,  showing  a  crisp  representation  of  the  discontinuity  with 
only  one  POD  mode  and  one  discontinuity  mode. 


2.2.2  Discontinuities  Capturing  using  Mathematical  Morphology 

The  POD  acceleration  methods  resulted  in  a  reduction  of  the  computational  time  of  two  orders  of 
magnitude  compared  to  the  full-order  model  (Pizmas  et  ah,  2008|),  as  mentioned  in  the  previous 
section.  In  spite  of  the  speed-up  values,  the  effectiveness  of  the  proper  orthogonal  decomposition 
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Figure  2.5:  Velocity  profile  at  t  =  1.0  for  full-order  model,  reduced-order  model 
and  no  discontinuity  mode,  and  reduced-order  model  with  2  modes,  including 
mode  (Premier  et  al.,  2009|). 


with  20  modes 
a  discontinuity 


method  is  challenged  by  the  presence  of  flow  discontinuities,  such  as  shock  waves  QLucia  et  al.|,  |2003[). 
To  improve  POD  robustness  and  performance  when  dealing  with  discontinuities,  such  as  shock 
waves,  we  developed  a  POD  method  in  which  the  basis  functions  were  augmented  by  discontinuity 
basis  functions,  as  shown  in  2.2.1.  The  discontinuity  basis  functions  were  generated  by  detecting 


the  flow  discontinuities  in  the  full-order  model  using  a  morphology  method.  This  section  presents 
the  methodology  we  used  to  detect  these  discontinuities. 

Developed  in  the  1960’s  by  Georges  Matheron  and  Jean  Serra  at  the  Ecole  des  Mines  de  Paris, 
France,  mathematical  morphology  is  a  method  for  the  analysis  of  spatial  structures  (|5oille|,  2003). 
Morphology  is  heavily  based  on  set  theory,  integral  geometry,  and  lattice  algebra.  Morphology  was 
originally  developed  to  describe  porous  media,  where  the  porous  material  can  be  characterized  as  a 
binary  image  containing  pores  or  solid  surrounding  material  (|Soillc ,  2003).  While  this  method  has 
primarily  been  used  in  the  image  processing  community,  it  also  has  been  used  on  a  wide  variety 
of  spatial  structures.  Herein  the  theory  and  application  of  morphology  is  extended  to  locating 
discontinuities  in  flows. 

Morphology  attempts  to  extract  structures  of  interest  located  within  the  whole  spatial  domain, 
herein  referred  to  as  an  image,  by  performing  various  set  operations  over  so-called  structuring 
elements.  Structuring  elements  are  used  to  probe  the  image  with  various  shapes  that  describe  the 
structure  of  interest.  Structuring  elements  can  take  on  any  number  of  shapes,  from  simple  rods 
(straight  elements)  to  circles,  or  in  extreme  cases,  galaxy  shapes.  Structuring  elements  can  also 
be  applied  in  combinations  to  create  larger  structuring  elements  of  varying  shapes.  While  there 
is  not  a  limitation  to  what  shapes  or  size  the  structuring  elements  can  take,  these  structuring 
elements  must  be  aligned  in  the  same  direction  throughout  the  image.  Failure  to  maintain  the 
same  alignment  throughout  will  distort  the  structures  that  are  sought.  In  addition,  structuring 
elements  can  only  be  applied  to  highly  structured  images.  Images  not  in  a  highly  structured  state 
must  be  transformed  to  this  state  before  applying  the  morphological  operations  as  usual. 
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The  structuring  element,  denoted  as  b,  is  a  locally  defined  quantity  in  the  local  coordinates 
( k,l ).  The  values  of  (k,l)  are  defined  depending  on  the  total  size  of  the  structuring  element.  If 
(0, 0)  is  the  point  of  interest  in  local  (k,  l)  coordinates,  then  for  a  five  point  vertical  rod  structuring 
element  the  global  coordinates  would  be  (i,j  —  2l),  ( i,j  —  l ),  (i,j),  ( i,j  +  l ),  and  (i,j  +  2l).  b  assumes 
the  value  of  the  image,  /,  at  these  coordinates. 

While  morphology  can  perform  a  wide  variety  of  operations  over  the  structuring  elements,  all 
of  these  operations  can  be  described  by  combinations  of  two  elementary  operations:  erosion  and 
dilation.  Both  of  these  operations  can  be  described  as  shifts  in  the  data  set.  Erosion  is  the  minimum 
of  these  shifts  of  the  image,  f(i,j),  over  the  structuring  element,  b(k,l),  given  as 

e(/)  =  min(/(i  +  k,j  +  1)  -  b(k,l)).  (2.19) 

Dilation  is  the  maximum  of  these  shifts  over  the  same  structuring  element, 


d(f)  =  max(/ (i  -k,j  -  l )  +  b(k,l)). 


(2.20) 


Using  equations  (|2.19|)  and  (2.20|)  in  conjunction,  more  powerful  operations  can  be  performed. 
While  these  operations  are  not  inherently  morphological  in  nature,  they  are  useful  for  detecting 
discontinuities  in  flows.  These  operations  include  blurring,  blur /erosion,  dilation/blur,  and  thresh¬ 
olding. 

Blurring  is  used  in  image  processing  as  a  noise  reduction  technique.  This  operation  is  required 
for  a  preliminary  elimination  of  extreme  values.  Blurring  is  also  required  for  locating  step-type 
edges  ([Lee  et  ah,  1987).  Without  blurring,  the  area  around  a  step  edge  is  not  a  local  maximum  or 
minimum,  and  the  erosion  and  dilation  operations  will  not  extract  discontinuities  of  that  type  (Lee 
1987).  Blurring  is  defined  as 


et  al 


ELD  m 

K 


blur(f)  = 


L  = 


K-  1 


(2.21) 


Here  blurring  is  defined  over  the  reduced  stencil  size  of  nodes,  K. 

While  blurring  is  necessary,  using  a  reduced  stencil  that  is  too  large  will  over-blur  the  image, 
resulting  in  overly  distorted  edges  that  are  impossible  to  discern.  A  good  reduced  stencil  size  for 
blurring  is  K  =  3  ([Lee  et  ah ,  1987);  i.e.,  when  using  a  vertical  rod  structuring  element  and  looking 


at  the  center  point,  the  values  used  to  blur  the  center  point  are,  in  global  coordinates:  (i,  j  —  l ),  (i,  j ), 
and  ( i,j  +  /).  A  consequence  of  this  stencil  size  for  blurring  is  that  the  structuring  elements  must 
be  at  least  two  elements  larger  than  the  blurring  stencil  in  order  to  mitigate  possible  over-distortion 
of  the  edges. 

The  blur /erosion  (be)  and  dilation/blur  (db)  operators  were  proposed  by  |Lee  et  al.  (1987|).  These 
operations  act  over  the  structuring  elements  and  are  used  to  eliminate  obvious  points  where  edges 
are  not  located  using  the  standard  erosion  and  dilation  operations  in  conjunction  with  blurring. 
These  combination  operators  initially  blur  the  image  and  then  apply  the  erosion  and  dilation 
operators  to  the  blurred  image.  The  be  and  db  operators  are  defined  as 


be(f)  =  blur(f)  -  e(blur(f )) 


and 


db(f)  =  d(blur(f ))  -  blur(f)., 


(2.22) 


(2.23) 


respectively. 
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The  last  operation  required  to  detect  edges  is  thresholding.  Thresholding  is  an  image  segmen¬ 
tation  technique  used  to  isolate  areas  of  interest  in  the  image.  For  non-binary  images,  thresholding 
is  almost  always  required  to  completely  extract  the  desired  information  from  the  image.  Grayscale 
images  can  take  on  a  wide  variety  of  values,  causing  the  previously  described  operations  to  assign 
more  than  two  points  to  an  edge,  where  the  edge  is  defined  as  being  only  two  points  in  width; 
thresholding  ensures  that  each  edge  is  only  two  points  in  width.  Thresholding  can  be  applied  glob¬ 
ally  or  semi-globally,  where  globally  is  a  binary  fix  to  the  image  and  semi-globally  maintains  the 
original  value  of  the  edge.  We  use  global  thresholding,  defined  as  (Ritter  and  Wilson|,  2000) 


thres(f)  = 


1  a<  f  <  P 

0  otherwise 


The  constants  a  and  /3  are  user  defined  values  specific  for  the  image  or  set  of  images  if  the  group 
of  images  describe  the  same  general  phenomenon. 

All  of  the  previously  described  operations,  when  combined,  allow  for  the  detection  of  edges  in 


images  with  two-point  edge  strength.  The  algorithm  for  edge  detection,  described  by  Lee  et  al. 


(|1987|),  is  adapted  here  for  multiphase  flow. 

1.  Read  in  the  images  with  discontinuities  to  be  detected. 

2.  Create  the  structuring  elements.  Herein,  we  used  four  rod-type  structuring  elements  that  are 
five  units  in  length  and  have  the  following  orientations:  horizontal,  vertical,  and  diagonal  (45 
and  135  degrees  with  respect  to  the  horizontal). 

3.  Blur  the  image  using  Equation  (2.21|)  with  K  =  3. 


4.  Compute  the  erosion  and  dilation  of  the  blurred  images  using  equations  ( 2.19)  and  (2.20). 

5.  Compute  values  of  the  blur/erosion  and  dilation/blur  operators  with  equations  (2.22)  and  ( 2.23) 

6.  Determine  the  maximum  values  of  the  be  and  db  operators 


ma x(bei(f)) 


i  =  1, ...,  4 


and 


ma x(dbi(f)) 


i  =  1, ...,  4. 


where  i  represents  each  of  the  four  structuring  elements. 

7.  Determine  the  edge  strength 

Indicator(f)  =  min(max(d6(/)),  max(6e(/))). 


(2.24) 


8.  Apply  the  global  thresholding  to  the  images  given  by  Equation  (2.24).  This  gives  edges  for 
the  discontinuity  that  are  only  two  points  in  width. 

The  method  developed  to  detect  discontinuities  was  tested  on  an  airfoil  operating  at  Mach  0.85. 
Figure  |2.6|  shows  the  edge  strength  index  (2.24)  contours.  This  indicator  captures  not  only  the 
shock  discontinuity  shown  in  Fig.  |2.7|  but  also  the  discontinuities  due  to  the  vortex  shedding  and 


the  interaction  between  the  shock  and  the  boundary  layer.  The  discontinuities  were  detected  with 
two-point  accuracy. 

Mathematical  morphology  was  shown  to  extract  key  features  of  the  flow  by  using  a  simple  set 
of  operations.  Regardless  of  the  original  grid  arrangement  of  the  data,  the  morphological  detection 
algorithm  captured  the  discontinuities. 
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Figure  2.6:  Edge  strength  index. 


Figure  2.7:  Mach  number  contours. 
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Chapter  3 

A  Novel  Nonlinear  Beam  Model 


This  chapter  presents  the  development  of  a  novel  nonlinear  beam  model.  This  work  was  motivated 
by  the  need  to  reduce  the  computational  time  for  structural  modeling  and  to  allow  us  to  use 
continuation  methods.  The  first  section  presents  the  derivation  of  the  nonlinear  equations  of  motion, 
followed  by  the  modal  representation  of  the  beam.  The  third  section  describes  the  numerical 
implementation  of  the  nonlinear  beam  solver,  followed  by  the  solver  validation.  The  last  section 
briefly  presents  the  development  of  a  finite  element  method  code  for  plates  that  was  needed  to  avoid 
using  canned  commercial  codes  and  to  maintain  control  of  the  source  code  while  doing  aeroelastic 
simulations. 


3.1  Derivation  of  the  Nonlinear  Equations  of  Motion 

3.1.1  Overview 

In  this  section,  the  nonlinear  flexural-flexural-torsional  equations  of  motion  are  derived  for  a  beam 
with  a  straight  elastic  axis,  along  which  the  cross  sections  can  vary  arbitrarily  and  abruptly.  The 
cross  sections  of  the  undeformed  beam  can  have  a  center  of  mass  offset  from  the  elastic  axis  and  the 
principal  centroidal  and  bending  axes  variably  and  uniquely  rotated  about  the  elastic  axis.  Finally, 
due  to  the  consistent  retention  of  nonlinear  terms,  the  equations  are  valid  for  large  deformations, 
within  the  scope  of  the  validity  of  elasticity. 


3.1.2  Definition  of  Parameters 


The  beam  is  assumed  to  have  an  elastic  axis  that  is  straight  when  undeformed.  The  inertial  reference 
coordinate  system  M  (x,  y ,  z)  is  located  at  the  fixed  end  of  the  beam,  where  the  x-axis  is  coincident 
with  the  undeformed  elastic  axis  and  positive  in  the  direction  of  the  free  end.  The  cross-sectional 
coordinate  system  B  (£,??,£)  is  fixed  to  each  cross  section  during  deformation.  When  the  beam  is 
undeformed,  the  cross-sectional  coordinate  system  is  parallel  to  the  reference  coordinate  system.  To 
describe  the  orientation  of  the  cross-sectional  coordinate  system  relative  to  the  reference  coordinate 
system,  three  Euler  angle  rotations  are  invoked.  The  inertial  frame  is  rotated  an  angle  i/j  about  the 
£-axis,  9  about  the  first  intermediate  y-axis,  and  cj)  about  the  second  intermediate  x-axis,  which  is 
in  turn  the  £-axis.  Each  of  these  angles  is  a  function  of  time  t  and  position  s  along  the  deformed 


elastic  axis.  These  coordinate  systems  and  transformations  are  similar  to  those  used  in  Crespo  da 
Silva  and  Glynn  (|1978a).  Figure  [hi]  illustrates  the  two  coordinate  systems. 
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Figure  3.1:  Diagram  of  Deformed  Beam 


To  characterize  the  beam,  the  following  parameters  are  used. 


m(s) 

=  j P  dp  d( 

(3.1a) 

jds) 

=  J f^P  {v2  +  (2)  dp d( 

(3.1b) 

Ms ) 

=  gpfdpdC 

(3.1c) 

j'c(s) 

=  JfAPV2dpdC 

(3.  Id) 

hds) 

=  -  JJ  pp(dpd( 

(3-le) 

D&) 

=  GK 

(3. If ) 

Dv(s) 

=  JJa  E  C2  dp  d( 

(3-lg) 

Dd*) 

=  J  j  Ep2  dp  d( 

(3-lh) 

Dnds) 

=  —  J j  Ep(dpd( 

(3.  li) 

In  this  report,  nonlinearities  will  be  described  by  an  order  that  will  indicate  how  many  factors  of 
the  displacements  u,  v,  w,  and  4>  comprise  a  term.  These  factors  may  include  spatial  and/or  time 
derivatives  of  the  displacements. 
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3.1.3  The  Lagrangian 

The  angular  velocity  expressed  in  the  cross-sectional  frame  is 


u(s,t)  = 


cp  —  ip  sin  6  'j  ( 

ip  cos  9  sin  i p  +  9  cos  <p  /  =  l  cov 
ip>  cos  9  cos  <p  —  9  sin  <p)  )B  v 


(3.2) 


Love’s  kinetic  analogy  flCrespo  da  Silva  and  Glynn ,  1978a;  Love|,  1944)  enables  the  curvature 
vector  to  be  calculated  in  a  manner  similar  to  that  of  the  angular  velocity,  replacing  temporal 
differentiation  with  spatial  differentiation: 


cp'  —  ip'  sin  9 

p(s,t )  =  ^  ip'  cos  9  sin  <p  +  9'  cos  <p) 
ip'  cos  9  cos  (p  —  9'  sin  cp) 


(3.3) 


One  of  the  fundamental  assumptions  of  the  beam  model  is  that  cross  sections  perpendicular  to 
the  elastic  axis  remain  perpendicular  and  rigid.  The  inertial  velocity  of  each  point  on  a  cross  section 
of  the  beam  can  be  decomposed  into  the  sum  of  the  velocity  of  the  intersection  of  the  elastic  axis 
with  that  cross  section  and  the  velocity  of  the  point  on  the  cross  section,  relative  to  the  intersection: 


1 

r  *  l 

i  i 

r  °  ] 

H 

* 

>  +  U)  X  < 

V  > 

i 

l w  J 

l  C  Ji 

Resolving  the  second  velocity  contribution  into  the  inertial  frame  provides 


V  = 


rjiv^)  cipc9  +  rjLU£  {scpsip  +  ccpcips9 )  +  (ccpsip  —  cipscps9 )  +  u 

r]LU()  sipc9  —  rjLu^  (scpcip  —  ccpsips9 )  —  (c cpcip  +  sipscps9 )  +  v 

—  ((u>ri  —  rjuiQ)  s 9  +  rju^ccpc6  —  (u^scpcd  +  w 


Correspondingly,  the  cross-sectional  kinetic  energy  is  determined  by 


T=  \  jj^  p  V  ■  V  dq  d( . 
which,  when  expanded,  simplifies  to 

T  =  ^ m  ( u 2  +  v2  +  w2)  +  ^  (jscu|  +  jv +  jvc^C 

+  me^uj^  [(c^s^  —  scpcips9 )  u  —  ( scpsips9  +  ccpcip)  v  —  scpc9w\ 
+  me^Lov  [cipc9u  +  c9sipv  —  s 9w\ 

+  merjUi^  [(s^sV’  +  ccpcips9)  ii  +  (ccpsipsO  —  scpcip )  v  +  ccpc9w\ 
—  mevLU ^  [cipc9u  +  c9sipv  —  s 9w\  . 


(3.4) 


(3.5) 


(3.6) 


(3.7) 


Additionally,  the  potential  energy  of  each  cross  section  is 

V  =  -  +  DVP^  +  ^c^c)  Dv(Pr)P(-  (3.8) 

With  the  kinetic  and  potential  energy  for  each  cross  section  known,  the  cross-sectional  Lagrangian 
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is  simply 


l  =  T  —  V.  (3.9) 

3.1.4  Inextensionality  Constraint 

It  is  assumed  that  the  length  of  the  beam,  specifically  the  elastic  axis,  remains  constant  during 
deformation.  Applying  the  Pythagorean  theorem  to  an  infinitesimal  length  ds  of  the  beam  yields 

ds2  =  dv2  +  dw2  +  (ds  +  du)2  (3.10) 

or 

1  =  v'2  +  w'2  +  (l  +  u')2  .  (3.11) 

Furthermore,  this  constraint  enables  i/j  and  9  to  be  expressed  in  terms  of  the  derivatives  of  the 
displacement: 


ijj  =  tan 


-l 


1  +  u' 


9  =  tan 


-i 


— w 


\J(l  +  u')2  +  v'2 


(3.12) 

(3.13) 


3.1.5  Hamilton’s  Principle 

With  the  Lagrangian  known  for  each  cross  section  and  multiplying  the  inextensionality  constraint  by 
a  Lagrange  multiplier,  the  extended  form  of  Hamilton’s  principle  provides  the  governing  equations 
of  motion  associated  with  the  variations  5u,  5v ,  Sw,  and  5cj)  (|Crespo  da  Silva  and  Glynn,  1978a; 
Meirovitch ,  1967]) .  These  are  derived  in  Appendix  |A|  and  presented  below: 


G'u  = 
G'v  = 
G'w  = 

Aa,  = 


A, 


d'tp 

du' 


89 


+  +  A  +  “0 


du' 

89 


d2l 


dtdu 


.  -  Qt 


.dip  w  , 

A*w  +  Aew  +  Xv 


d2l 


dtdv 


Qv 


'  89  .  ; 

Ae^—,  +  Aw 
dw’ 


d2l 


where 


,  an  an  ai  ,  ,, 


(3.14) 

(3.15) 

(3.16) 

(3.17) 

(3.18) 


The  above  equations  are  simplified  by  carrying  out  the  following  substitutions  and  Taylor  series 
expansions: 


if)  =  tan 


-l 


/  /  /2  /2  ' 
V  ,  V  W 

( 1  +  IT  +  T" 


(3.19a) 
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dij} 

du' 

dip 

dv' 

6 

d6_ 

du' 

dv' 

d6 

dw' 


(1  +  u')2  +  v'2 

(1  +  0  _  v'2  w' 

(1  +  u')2  +  v'2  ~  T  +  ^ 


(3.19b) 

(3.19c) 

(3.19d) 

(3.19e) 

(3.19f) 

(3.19g) 


u 


/ 


u  = 


(3.19h) 

(3.19i) 


The  terms  in  the  extensional  equation  (3.14)  are  expanded  in  Taylor  series  to  quadratic  order,  and 
the  terms  in  the  bending  and  torsional  equations  Q3.15)-(3.17)  and  expanded  to  cubic  order.  The 
extensional  equation  only  requires  quadratic  order  as  it  is  used  to  determine  A,  which  only  appears 
when  multiplied  by  a  spatial  derivative  of  a  displacement  in  the  bending  equations. 


3.1.6  Application  of  the  Galerkin  Method 

Invoking  the  Galerkin  method  facilitates  discretization  of  the  beam,  particularly  with  regard  to 
abrupt  property  changes  along  the  length.  For  a  linear  cantilevered  beam,  the  boundary  conditions 
are  known  (pisplinghoff  et  a! ,  |1996).  For  the  case  of  the  nonlinear  beam,  the  boundary  conditions 
remain  the  same  and  are  obtained  from  the  boundary  components  resulting  from  integration  by 
parts  of  Hamilton’s  principle  (Crespo  da  Silva  and  Glynn,  1978b). 


3. 1.6.1  Bending  Equations  of  Motion 


The  following  set  of  mode  shapes  QBisplinghoff  et  ah,  1996)  is  used  when  applying  the  Galerkin 
method  to  the  bending  equations  of  motion. 


IF,;  (s)  =  cosh 


where 

cosh  fli  +  cos  Pi 

<7j  =  - 

sinh  pi  +  sin  Pi 
1  +  cos  Pi  cosh  Pi  =  0. 


Additionally,  Wi(0)  =  W/( 0)  =  W('{L)  =  W["{L)  =  0. 


(3.20) 


(3.21) 

(3.22) 
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Multiplying  ([3. 15)  and  (3.16|)  by  W,;.  expanding  the  rightmost  side,  integrating  over  the  length 
of  the  beam,  and  performing  integration  by  parts  yields 


Wi(L) 


.dil)  A  d6 

A*avi+A‘dv' 


f 


S  =  L 


=  /  Wi 

/  o 


mv  —  me„  (  <p> +  (pep  +  cpv  w  +  2cpw  v  +  v  +  2 cpv  w  +  2(pv w‘ 


—  mejj  (y'v'  +  cpw'v'  +  (pv'w1)  +  me(  ^ '  cp 4>2  —  (p  +  +  2(pv'v'\ 


+  me^  ( cpv '2  —  2 v'w'  +  cpv'v'  —  w'v'  —  v'w')  —  Qv 


ds 


and 


Wi(L) 


A 


dO 

dw' 


=  [  Wi 
Jo 


-  s=L 


mw  —  me^ 


p2  —  <p  +  ^<P2&  +  ^4> w ,2  +  2  cpw'w' 


—  me, 


(i pw'2  +  cpw'w')  —  me £  (j)2  +  cp <j>  +  w 12  +  w'w'^j  —  Qv 


ds. 


(3.23) 


(3.24) 


3. 1.6. 2  Torsional  Equation  of  Motion 


The  set  of  mode  shapes  ( Bisplinghoff  et  al.,  1996)  used  when  applying  the  Galerkin  method  to  the 
torsional  equation  of  motion  is 


$i(s)  =  \/~2  sin  ^Jj- , 


where 


7i  = 


(2i  —  1)  7T 


(3.25) 


(3.26) 


Furthermore,  4>j(0)  =  4>'(L)  =  0. 

Analogous  to  the  bending  equations,  (|3.17)  is  multiplied  by  and  integrated  over  the  length 
of  the  beam,  providing 


[  A^ids  =  [  Q^ids. 
/  0  Jo 


(3.27) 


3.1.7  The  Lagrange  Multiplier 


The 


terms  of  interest  in  ([3.23|)  and  (3.24)  are  /  Wi  (A?/)  ds  and  /  Wi  (A w')  ds ,  which 

Jo  Jo 

p— .  a  dtp  .  90 


are 


obtained  from  (3.14).  Letting  Au(s,t )  =  A^,——  +  Ag  — — -  and  observing  that  u  and  the  time  and 
-  dw  dw 
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space  derivatives  thereof  are  quadratic  terms,  a  quadratic  expression  for  Au  is  obtained: 

Au(s,  t )  =  (v'w'  +  v'w')  -  j^v'v'  -  jvw'w ' 

+  D't-v'v"  +  D'^w'w"  +  D^v'v'"  +  D^w'w'" 

—  D'^  (w'v"  +  v'w"')  —  Dr^  (w'v'"  +  v'w'")  . 

Equation  (|3. 14 )  is  integrated  from  the  free  end  of  the  beam  to  s: 


(3.28) 


/; 


[AM(s, t)  +  A(s, i)  (l  +  t/)]7 ds  =  j  mu  +  mev  <j>w'  —  2<j)w'  —  v'  —  fiw'^J 

+  me(  f  cfyv '  +  2 <jyb'  +  cj>v'  —  w'j  —  Qv 


ds 


(3.29) 


Since  A  is  quadratic,  (1  +  u')  1  «  1.  Knowing  this,  that  Gu(L,t )  =  0  (Crespo  da  Silva  and  Glyim|, 
1978b),  and  invoking  the  expression  for  u  obtained  in  (3.19i[),  (|3.29|)  simplifies  to 


—  m 


A  (s,t)  =  - Au(s,t )  + 


+rae(  \  4>v'  +  2 <j)v'  +  cj>v'  —  w'j  —  Q 

and  at  L, 

A  (L,  t )  =  [-jtf  (v'w'  +  v'w')  +  jcy'v'  +  jvw'w']  s=L  . 


f  ( i )'2  +  v'v  +  w'2  +  w'w)  ds 

Jo 


+mev  (  —cj)w'  —  2 cj)w'  —  v'  —  (j)w 


ds, 


(3.30) 


(3.31) 


rL  rL 

,  /  Wi  (At/)  ds  and  /  W,;  (A w')  ds  are  determinable.  Beginning  with  integration 

Jo  Jo 


r  L  rL 

With  A  known,  /  Wi  (A?/)  ds  and 

'  o  .70 

by  parts  of  the  former, 

r -L  rL 


[  Wi  (Xv'Y  ds  =  —  [  W'v'\ds  +  Wi(L,t)\(L,t)v'(L,t), 
Jo  Jo 


(3.32) 


which  is  expanded  as 

rL  rL 


f  Wi  (Xv1)' ds  =  f  W'v' Au(s,t)ds 

Jo  Jo 

—  f  W'v1  f  <  —  m  f  (v'2  +  v'v  +  w'2  +  w'w')  ds 

Jo  Jl  {  Jo 

+  mev  ^ — (f)w'  —  2 <fiw'  —  v'  —  (f)w''sj 

+  meQ  (j>v'  +  2 (j>v'  +  cj>v'  —  w'^j  —  j-dsds 
+  Wi(L,t)X(L,t)v'(L,t). 


(3.33) 
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Au  contains  derivatives  of  the  spatially  dependent  parameters,  which  are  eliminated  by  further 
integrating  by  parts,  ultimately  yielding 

^  Wi  (A v')'ds  =  -  f  |dc  (v'v"2Wl  +  v'2v"'W'  +  v'2v"W"  +  v'v"2W^ 

+  Dv  (v'w"2W[  +  v'w'w"W[  +  v'w'w"W"  +  v"w'w"W[ 

-  DnC  ( v"2w'W'  +  v'v"'w'W'  +  v'v'w'W" 


-  D,nC  \3v'v"w"W(  +  v'2w"'W[  +  v'2w"W”j  jds 

+  J  [v'w'  +  v'w')  —  j^v'v1  —  j'qw'w' 

+  Dqv'v'"  +  Drjw'w"1  —  Dv£  [w'v"'  +  v'w'")  |ds 

—  f  W[v'  f  |  —  m  f  [y'2  +  v'v'  +  w'2  +  w'w')  ds 

Jo  Jl  {  Jo 


+  mev  y—tfiw'  —  2 <fiw'  —  v  —  4>w' 

+  me(  ( <j)v'  +  2(f>v'  +  (f>i )'  —  wj  —  Qu  >dsds 


+  Wi 


./  ,  „./2  ../ 


/2  ■•/ 


JrjQ  I  V  V  W  +  V  W  I  +  J(V  "V  +  JVV  W  W 


-  S  =  L 


and  similarly 
r -L 


I 


Wi  (A  w')'  ds  =  -  j  |d?  (w'v"2W[  +  w'v'v'"W[  +  w'v'v'W"  +  v'v"w"W^ 

+  Dv  (w'w"2W[  +  w'2w'"W[  +  w,2w"W"  +  w'w"2W'^j 

—  ^  v'  w"2W[  +  v'w'w"'W[  +  v'w'w"W"^j 

—  DvC  (liv"w'w"W'l  +  v"'w'2W[  +  u/Vtff  )  jds 

+  J  W-  w'  |  ( v'w '  +  v'w ')  —  jyv'v'  —  jvw'w' 

+  Dqv'v"'  +  Dvw'w'"  —  D ^  [w'v'"  +  v'w'") 

—  f  W[w'  f  |  —  m  I  ( i b'2  +  v'v'  +  w'2  +  w'w')  ds 

Jo  Jl  l  Jo 

+  mev  (  —  j>w'  —  2(j)w'  —  v'  —  (f>w' 


+  me(  ( <pv'  +  2 +  (f>v'  —  w')  —  Qu  >dsds 


+  Wi 


,/2  .  /  /  ../ 


/2 


—jri(  [v  w  +  v  w  iv  j  +  j^V  V  w  +  jvw  U! 


(3.34) 


-  s=L 


(3.35) 
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3.1.8  The  Equations  of  Motion 

At  this  point,  everything  can  be  determined.  Remaining  terms  containing  derivatives  of  the  spatially 
dependent  parameters  are  integrated  by  parts  to  distribute  the  derivative  into  the  co-factors. 
Consequently,  the  equation  of  motion  in  the  y-direction  for  each  i  is 


Wi(L)  (me,  —  /  {v'2  +  v'v'  +  w'2  +  w'w1')  ds  +  me,  (— vv'  —  cj)vw'^ 

L  Jo 

+me(  ( (j)vv‘ '  —  vw')  +  jg  (^< f>w 1  +  <j)w '  +  2 w'v'w'  +  w'2v'^j  +  jv  i^—(j)w'  +  cj)2v' 
+jv  (y-(j)w'  +  2 4><j)V,'SJ  +  ^  —  2(f)(j)v'  +  v'v'2  +  cj>w'  +  v'w12  +  v'  —  (j)2v' 

+j)Q  ^ v'2v '  +  4>w'  +  v'w'w'^j  +  jvQ  ( 2(j)v '  +  Acjxfrw1  +  2(j)v'  —  w'  +  2 cf)2w' 

—  -jv(  (VV  +  u/V)  —  D^w'cj)1 


s=L 


W-  ( mev  —  me^(j))  /  ( v '2  +  v'v'  +  w'2  +  w'w')  ds 

l  Jo 

-\ -me,  {—vv'  —  (pvw ')  +  me^  {4>vv'  —  vw')  +  j £  (q iw '  +  <j>w'  +  2 w'v’w'  +  w'2v' 
+jr?  iy-^w'  +  4>2v  —  4>w'  +  2 (jxj)v''\  +  jq  ^—2 <j><j)v'  +  v'v'2  +  cjytb'  +  v'w'2  +  v' 
+JC  +  v'2v'  +  +  v'w'iv'^j  +  jvq  (2 <j)v'  +  4 <j)<f)w'  +  2cj)v'  —  w' 

1  1/2  -.A  (  j  j/2  .  //„/.// 


+jri(  I  2(j)2w'  —  -v'*w'  —  -w'*w'  )  +  Dc  (  v'v"“  +  v"w'w 
+Dvq  (^—v'v"w"  +  w'w"2^j  +  v'  J  —m  J  {v'2  +  v'v1  +  w'2  +  w'w )  ds 
+me,  (—(f>w'  —  2(j>w  —  v'  —  (ftw'^J  +  meq  (jf>v'  +  2 <pv'  +  (jiv  —  w'^j  —  Quds 
D ^  (yj)'w'  +  w'2v"^j  +  {( f>2v' '  —  cf>w")  +  (y"  —  4>2v"  +  v'2v"^j 

[(pw"  +  v'w'w"')  +  Dvq  ^2cj)v"  —  w"  +  2 cf)2w"  —  \w'2w"  —  ^w'2w'^j  |ds 
J  Wi  (ynv  —  mev  (j)2  +  4><j)  +  <j>v'w'  +  2 fyw'v’  +  v'2  +  2 cj)v'w'  +  2 4>v'w 
—merj  {+v'v'  +  <j>w'v'  +  <j)v'w ')  - 

+me^  (2 cj)v'v'  +  cl)v'2  —  2 v'w'  +  <j>v'v'  —  w'v'  —  v'w' )  —  Qv  )  ds. 


+W'' 


b2  -  J>+  \<t>V2 


(3.36) 
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Additionally,  the  equation  of  motion  in  the  ^-direction  is 


Wi(L) 


(■ me^cj)  +  me()  /  ( v '2  +  v'v'  +  w  2  +  w'w')  ds  —  mev  ( cf>vv ’  +  4>ww') 

Jo 

+me(  (— ih/  —  ww')  —  (<fiv'  +  w'v'2^)  +  j v  i^—j>v'  —  2 <p<j)w'  —  <\>v'  +  w'^J 
+jv  ^ +  w'w 12  +  w'2w'^j  +  j ^  (j)v'  +  w'v'2  +  2 (fxfiw1  +  4>v'  +  (f)2w ' 
+j^v'w'v'  —  jv£  (—4 4><])v'  +  v'v'2  +  2 4>w'  +  2 w'v'w'  +  t/ii/2  +  i/  —  2<j!)2'i3/ 


•  /  /2  ../  .  f  /2  ../  !  C\  1  "  /  i  r,  /  /  ■■ , 

—Jr](  (  TjV  V  +  -lO  V  +  Z(pW  +  Zv  W  W 


s=L 


W' 


(mev4>  +  me^)  /  (i>/2  +  +  w'2  +  w'w')  ds 

Jo 

-me,  (<; t>vv'  +  i cj)ww')  +  me(  (—i)v'  —  ww')  —  ypv1  +  w'v'2')  —  jv4>v' 


+jv  (  — 2<^ii/  —  4>v'  +  w'  —  4>2w'  +  w'2w'  +  w'w'2  )  +  ( 4>v'  +  w'v' 


-/2 


+k 


+  4>v'  +  (jrw  +  v'v 'w'j  —  jr!(  y—44>j>v'  +  v'v  2  +  2 <j)w' 
1  ,2..,  ,  1.../2,./ 


—jv(  I  2 w  v  w  +  v  w  +  v  —  2 <p  v  +  -v  v  +  -w  v  +  2 <j)iv  +  2v  w  w 
\  Zj  z 

i  j~~\  ( jl  n  i  /  //2\  /  //2  .  T— ,  /  //  //  i  \  //////,  /  //2\ 

( (p  v  +  w  v  I  +  LV/U;  w  +  D^v  v  w  —  I  v  w  w  +  v  w  I 

ps  ps 

+w'  /  —  m  /  ( i )'2  +  v'v'  +  w'2  +  w'iv')  ds  —  me^iftw' 

Jl  Jo 

+me7)  ^—2 cj)w'  —  v'  —  cj)w')  +  me^  (j. w '  +  2 <j)v'  +  4>v'  —  w')  —  Quds 
Dv  +  w"  —  cj)2w"  +  w,2w"^j  +  (cf>v"  +  cj)2w"  +  v'v"w') 


+Wi 


(  —v"  +  2 cj)2v"  —  -v'2v"  —  -w'2v"  —  2cf)w"  —  2 v'w'w' 


>ds 


Wi  ( mw  —  mer. 


-me, 


f 

(4>w'2  +  <t>w'i v')  —  me^  (j)2  +  4><j)  +  w'2  +  w'iv'^j  —  Qw ^  ds. 


4>2  -  <f>  +  +  \  j)W'2  +  ^j>w'w' 


(3.37) 
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Finally,  the  torsional  equation  of  motion  is 


<I>, 


( mevw '  —  me^v')  /  (y'2  +  v'v'  +  w'2  +  w'w')  ds 

Jo 


1 


1 


+me„  — cf>v  +  iv  —  -(/)  w  —  v v  w  —  -ww 


+me(  (  —v  +  —  ^  +  2^V’2 


+j£  +  v'w'  +  w'v'^j  +  jrj  {  —  +  v'w1  +  0w'2) 

+j(  {fa'2  —  4>w12  —  v'w  )  —  jv£  (v12  —  w'2  +  4 (frv'w1) 
+DV  (jiv"2  —  v'w"  —  4>w"2^  +  ^—c/)v"2  +  v"w"  +  4>w"2^j 

+Dr^  (v"2  +  4 4>v"w"  -  w"2^j  +  (yf)'  +  v"w')  |(is  =  j  Q^ids. 


(3.38) 


The  effects  of  including  all  the  cubic  non-linear  terms  was  assessed  on  the  Heavy  Goland  wing. 
The  parameters  of  the  Heavy  Goland  wing,  which  include  airfoil  shape,  mass  properties  and  wing 


stiffness  are  presented  in  (Eastep  and  Olsen,  1980)  and  will  not  be  repeated  in  here. 

An  analysis  of  the  Heavy  Goland  wing  using  a  nonlinear  structural  model  and  a  simple  aero¬ 
dynamics  approach  was  presented  in  (ptrganac  et  ah,  2005|).  The  same  nonlinear  structural  model 
was  later  coupled  with  a  RANS  model  (|Cizmas  et  ah,  2007|).  The  structural  model  included  the  me¬ 
chanics  of  both  the  traditional  out-of-plane  and  torsional  components  plus  an  additional  in-plane 
mode  of  motion.  This  additional  mode  and  the  presence  of  nonlinear  effects  led  to  remarkable 
nonlinear  dynamic  interactions  (|Cizmas  et  al.,  2007).  In  both  cases,  the  structural  model  did  not 
include  all  the  cubic  nonlinear  terms. 

The  results  shown  in  Figs.  3.  Ml  highlight  the  effect  of  the  cubic  nonlinear  terms  of  the  struc¬ 
tural  model  and  underscore  the  importance  of  the  nonlinear  effects  in  aeroelastic  simulations.  Fig¬ 
ure  T2  shows  that  missing  some  of  the  cubic  terms  did  not  affect  the  LCO  prediction.  Figure-|T3|, 
however,  shows  that  the  numerical  simulation  using  all  cubic  nonlinearities  predicted  flutter,  while 
the  case  that  was  missing  some  of  the  cubic  nonlinear  terms  did  not.  Consequently,  in  the  flutter 
case,  it  was  important  to  keep  all  the  cubic  terms. 


3.2  Modal  Representation 

3.2.1  Introduction 


The  equations  of  motion  (|3.36|)-(3.38D  are  nonlinear,  coupled  partial  differential  equations  in  time 
and  space.  As  mentioned  previously,  the  boundary  conditions  are  known.  Therefore,  it  is  possible 
to  represent  the  deformations  as  an  infinite  sum  of  the  product  of  spatially  dependent,  orthogonal 
mode  shapes  and  time-dependent  coefficients.  Taking  the  mode  shapes  invoked  in  the  Galerkin 
method,  the  displacements  can  be  expressed  as  follows: 


l  m  n 

w  =  m(t)Wi(s)  V  =  Vi(t)Wi(s )  4>  =  4>i(t)$i(s).  (3.39) 

i= 1  i= 1  2=1 

Upon  invoking  these  mode  shapes,  the  partial  differential  equations  reduce  to  a  system  of 
ordinary  differential  equations  in  time,  for  which  the  time  coefficients  are  the  unknown  quantities. 
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0.0  0.5  1.0  1.5  2.0  2.5 
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(a) 


(b) 


Figure  3.2:  Heavy  Goland  Wing  time  coefficients  at  LCO:  (a)  out-of-plane  bending,  and  (b)  pitch¬ 
ing. 


Therefore,  the  equations  of  motion  are  written  in  the  following  manner  using  linear  and  non-linear 
arrays: 


[Ml  +  Mnl\ 


+  [Cl  +  Cnl]  |  v  |  +  [Kl  +  Knl\ 


Fl  +  Fnl, 


(3.40) 


where  w,  v,  and  <f>  each  represent  vectors  containing  the  time  coefficients:  {w\, ...,  wi}T ,  {ui, ...,  vm}T , 
and  { <7^1 ,  ...,4>n}T- 

Each  of  the  matrices  is  (l  +  m  +  n)  x  (l  +  m  +  n)  and  can  be  expressed  as  the  concatenation  of 
9  submatrices  with  the  following  dimensions: 


[l  x  l]  [lx  m\  [l  x  n] 
[m  x  l]  [m  x  m\  [m  x  n] 
[n  x  l]  [n  x  m\  [n  x  n] 


(3.41) 
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Figure  3.3:  Heavy  Goland  Wing  time  coefficients  at  flutter:  (a)  out-of-plane  bending,  and  (b) 
pitching. 


3.2.2  Linear  Matrices 
3. 2. 2.1  Linear  Mass  Matrix 

The  linear  mass  matrix  is  defined  as  follows: 

[ML}id  =  [L  (WimWj  +  W'jvW ')  ds  - 
J  0 

[Mdif+j  =  -  [  W!j„cW’ds+  [Wd,(w;i_L 

J  0 

[ML\i+m+j  =  I  Wime^jds 

J  0 

[ML}l+iJ  =  -  [LW'jv(;W'ds+[Wdv<W'}s=L 
J  0 

Wl+i,I+j  =  f  (WtmWj  +  W'jcW')  ds  - 

J  0 

[M4+y+m+i  =  -  [  Wime^jds 
J  0 

[ML}i+m+iJ  =  I  $ imerjWjds 

J  0 

[ML]i+m+i  i+j  =  -  [  $ imecWjds 
Jo 

[ML}l+m+il+m+j  =  I  §d&jds. 

J  0 


(3.42) 
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3. 2. 2. 2  Linear  Damping  Matrix 

In  the  equations  of  motion,  there  are  no  linear  occurrences  of  the  first  time  derivative  of  any  of  the 
deformations.  Therefore,  the  linear  damping  matrix  is  zero: 

[Cl]  =  0.  (3.43) 

3. 2. 2. 3  Linear  Stiffness  Matrix 

The  nonzero  submatrices  of  the  linear  stiffness  matrix  are 
[Kr,h=  j^Wl'D^ds 
[KL\iM  =  -  fwl'D^W'jds 

J  0 

[KL\l+i.  =  -  W'lD^ds 

Jo 

[tfdi+U+j  =  / W?D(Wfd, 

J  0 

[KL]i+m+ij+m+j  =  [  KD&ds.  (3.44) 

•J  0 

3.2.3  Nonlinear  Matrices 

In  determining  the  nonlinear  matrices,  there  are  several  possibilities  concerning  where  to  assign 
terms  from  the  equations  of  motion.  This  decision  is  done  systematically  such  that  terms  containing 
second-order  time  derivatives  are  assigned  to  the  mass  matrix,  remaining  terms  that  contain  first- 
order  time  derivatives  are  assigned  to  the  damping  matrix,  and  terms  containing  no  time  derivatives 
are  assigned  to  the  stiffness  matrix.  Furthermore,  for  terms  containing  factors  that  have  the  same- 
order  time  derivative,  placement  within  the  matrix  is  done  such  that  the  term  with  the  lowest- 
order  spatial  derivative  is  the  unknown.  Finally,  for  instances  in  which  the  temporal  and  spatial 
derivatives  of  the  factors  are  both  of  the  same  order,  the  spatial  derivative  of  v  is  taken  to  be  the 
unknown. 
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3. 2. 3.1  Nonlinear  Mass  Matrix 

The  nonlinear  mass  matrix  is 

cL 


[Mj\ n]ij  =  j  jw/  (me^cp  +  me^)  J  w'W'jds  +  W[  [— me^cpw'  -  me^w1]  Wj 
+  W'w '  J  —  rae()  Wjds 

+  W'  (-2jvc(p  -  2jvCv'w'  -  j^2  +  jvw'2  +  j(p2^j  W' 

—  W-w'  J  m  J  w'Wjdsds  +  Wi  [—me^pw1  —  me^w']  Wj  |ds 
—  Wi  (mevp  +  me()  /  w'Wjds  +  Wi  [—me^cpw1  —  me^u/]  Wj 

J  o 

+  Wj  (-2^0  -  2 +  jr,w'2  -  J77<^2  +  Jc^2) 


-  s=L 


[MNL\ij+l  =  J  |  W[  ( mev(j)  +  me^)  J  v'W-ds  +  W/  [-me^pv1  —  me^v']  Wj 
+  W'w1  J  (-me,  +  me^p)  Wjds 

+  W/  -  \hCv'2  ~  \hc,w'2  ~  +  k &  +  wj 

—  W-w'  J  m  J  v'Wjdsds^ds 

—  Wj  (mer/p  +  me()  /  v'Wjds  +  W*  [— mevpv'  —  Wj 

J  o 


+  w*  (  2 j^</>2  -  2.?^'(/2  -  2:jri(W,z  -  jvP  +  j'c0  +  kv'w' )  w: 


,./2 


-  s=L 


[Matl] 


j,Z+m+j 


W/u/  J  +  me^v')  &jds 


+  Wi 


1  ,o  1 


—me^  I  p-p>z  +  p- w ,2  )  -  me^ 


Tj  >  (is 
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iMNL\i+i  j 


[MNL\l+ij+j 


[Mnl] 


[MNL]l+m+iJ 


[Mnl] 


W[  (me,  —  /  w'Wjds 

Jo 

+  W[v'  J  {—men(j)  —  me()  W'jds 

rs  rs 

+  Wi  \—mev(f)v'  —  me^t/]  Wj  —  W[v'  m  w'Wjdsds 

Jl  Jo 

+  w'  ~  \  jv(v'2  ~  \jv(w'2  ~  Jv $  +  J +  kv'w>')  Wj  jds 

Wi  (mev  —  me^4>)  /  w'Wjds  +  Wj  [— mev(j)v'  —  me(U/]  W ) 

Jo 


+  Wi  (  2 j^c/)2  -  72:ir,Qv'2  -  ^jr,Cw>jl  ~  J +  J +  J(V'W'  )  Wj 


a 


-  s=L 


W-  (me,  —  mec4>)  /  v'Wjds 

10  l  Jo 

+  W/  [— (t/  +  cjnv')  +  me(  [cj)v'  —  u/)]  Wj 
+  WjV  J  (—mev  +  me^(j))Wjds 
+  Wj  [— me,,  (u'  +  </>u/)  +  me^  ((ftv'  —  w')]  Wj 

+  (2.7^  +  j?^/2  +  j n4>2  -  iU2  +  icu/2)  Wj 

—  W'v'  J  m  J  v'Wjdsds^ds  — 

+  Wj  [— mev  (y'  +  (j>w')  +  me^  (tyv'  —  w'  )]Wi 


Wj  (me,  —  me^4>)  /  v'W'As 
Jo 

—  w')  ]  Wj 

+  Wj  (2jv^  +  jyw'2  +  jv4>2  -  jc4>2  +  j(V,2^j  Wj 


j  |  WjV  J  (y—me^w'  +  me^v')  &jds  +  W'jyivUj 


—mev  (cj)  +  v'w ')  +  me^  (  -cj)2  +  -v'2 


+  w 

-  [W Jtw%\s=L 


}ds 


J  $i|  ( me^w'  —  me^j  J  w'Wjds 


+ 


—mev  I  -d>2  H — w'2  I  —  me. 


Wj  Us 


j  <3?i|  (me,w'  —  me^v')  J  v'Wjds  +  jyw'W'j 


+ 


—mev  (cj)  +  v'w')  +  me(j  (  -(j)2  +  -v'2 


Wj  Us. 


(3.45) 
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3. 2. 3. 2  Nonlinear  Damping  Matrix 

The  nonlinear  damping  matrix  is 


[CNL\ij 


[CNL]iti+j 


[CNL\itl+m+j 


[CNL\i+ij 


W[  (mev<f>  +  me()  /  w'Wjds  +  Wi  (—me^w'  —  me^ui')  Wj 
Jo 

+  W[  (—JnCuw  +  jrj'w'w  )  Wj  —  W[w  j  m  j  w'Wjdsds | ds 
Wi  (me^  +  me^)  f  w'Wjds  +  Wi  {—j^v'w'  +  jvw'w')  Wj 

Jo 


W[  ( men<j>  +  me^)  J  v'Wjds 
+  W[  {—jnCv'v'  ~  ^Jricw'w'  —  jgw'v'  +  j(v'w')  Wj 

—  W-w'J  m  J  v'Wjdsds^jds 

—  Wi  (merjip  +  me()  /  v'Wjds 

Jo 

+  Wi  {-j^v'v'  -  2 j^w'w'  -  j^w'v'  +  j(v'w')  Wj 


s=L 


io 

+  Wn 


W\w  J  (—2 mevw'  +  2 me^v')  &jds 
-me,  {^4>(b  +  2 w'w'^j  —  me^cj)  &j  +  W'  —  2jv^w')  &j 

+  wi  -  Jv'ij>  -  ^J-n^w  +  jcy’  +  2 jrfw')  ®j  jds 

-  [Wi  (4 jnrfv'  -  2 jvCw'  -  jty'  -  jvv'  -  2 jv<jyw'  +  j(v)  <t>,]  g=L 
-2  [Wij^w'^.j\s=L 


W[  (me,  -  metf)  [  w'W'ds  +  W^j(v'w'W' 

Jo 

W[v'  J  m  J  w'Wjdsds^jds 

Wi  ( me,  —  me^(f>)  f  w'Wjds  +  Wijcy'  w'W'j 
Jo 


-  s=L 
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[CNL\i+iti+j  =  |  w'  (mev  -  me((j))  J  v'W'ds 

+  Wi  [—mev  (y1  +  2(/>u/)  +  me^  ( cj)v '  —  2 w/)]  Wj 
+  W[  {2j^w'w'  +  j^v'v')  Wj  —  W[v'  J  m  J  v'Wjdsds^ds 

—  Wi  (me,  —  me^cj))  f  v'Wjds  +  Wi  (2 j^w'w'  +  j(v'v')  Wj 

Jo 


s=L 


[CNL\i+ij+rn+j  =  J  |  Wy  (—2meriwl  +  2mec«')  $jds 

+  Wi  -me,  (j)  +  2 w'v'  +  2 v'w'^j  +  me^  (jicj)  +  2v'v'  ’ 

+  wi  {ZjriC.v'  +  4 j^w'  +  jtrw'  -  jvw'  +  2 jn(j>v'  +  jcu/)  <h, 

-  2W/^#'<h,- jds  -  [Wi  (2jri^v  +  4jnC<l)w  +  j^w  -  jvw ')  ®j\s=L 
~  [Wi  (2 jv(j)v'  +  j^w'  -  2 jC(jm')  ®j\s=L 


[ Cnl ] 


&i  ( mevw '  —  me^v')  J  w'Wjds 


+  4>i  {j^w'  +  -  j(4>w)  Wj  ^ds 


[Cnl] 


4>i  ( menw '  —  me^v')  J  v'W'jds  +  4>i  {—jr^v  ~  ^irt^w')  Wj 
+  4>i  (jgw'  -  jv4>v'  +  jrfjj  +  j(;v'<f)  -  Jiqw')  W'Xds. 


3. 2. 3. 3  Nonlinear  Stiffness  Matrix 

The  nonlinear  stiffness  matrix  is 

[KNL\ij  =  £  |  Wi  (DiV"2  +  Dvw"2  -  D^v"w")  W'j 
+  W['  £vw'w"  -  ^Dv(v"w'^J  Wj | ds 

[KNL\l+j  =  £  |  W[  (. Dft/V '  -  DvCw"2)  W'j 

+  W"  (d(v"w'  -  ^D^v'v"  -  2Dv(w'w"j  Wj  Ids 


(3.46) 
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[KNL]<l+m+j  =  £  [w'D,v"*i  +  W”  (-£>„«"  -  D^w"  +  Dcv"  +  D((j>w")  4>, 

+  W"  (2 Drjrfv"  -  2 D^w")  <S>£ds 

[Knl]1+1,  =  £  {  W[  ( Dcv"w "  +  DvCw"2)  Wj 

+  W. ■  [d^w'v"  -  ^ DvCw'w ")  W'£s 

[KNL]l+i,l+j  =  £  { Wi  (Dcv"2  -  D^’w")  Wj 

+  W. ■  [d^v'v"  +  D(w'w"  -  J),lCr'u£  Wj\ds 

[KNL]l+i>l+m+j  =  £  |  W'l  {D.n(j)v"  -  Dvw"  -  Drfv"  +  Dcw"  +  2 D^v")  <t>5 

+  2 Wl'Drftfni/'Qj  +  Wl'D^w'&^ds  +  [WiD^w'^'}  s=L 

[KNL]l+m+iJ  =  £  {-^D,qw"W;  +  Q'Ds v"W '}  ds 
J  0 

[KNL]l+m+i,l+j  =  £  ^  ( ~D,w "  +  Dcw"  +  D^")  Wjds 

J  0 

[KNL\l+m+iHm+j  =  £  ^  (Dy2  -  Dvw"2  -  D(y"2  +  Dcw"2  +  4 D^v'W')  ds .  (3.47) 


3.2.4  Forcing  Vectors 
3.2.4. 1  Linear  Forcing  Vector 

The  linear  forcing  vector  is 

—  /  W^iQwds 


{ Fl },  =  /’ 
Jo 

nyi+i  =  [ 

Jo 


WiQvds 


{Fl} ’l+m+i  ~  /  ^iQ(j)ds. 

Jo 


(3.48) 
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3. 2. 4. 2  Nonlinear  Forcing  Vector 

The  nonlinear  forcing  vector  is 


{FNL\i  = 
{Fnl\i+1  = 


[L  WiW'  f1 
Jo  Js 

[L  Wiv'  [L 
Jo  Js 


Qudsds 


Qudsds. 


(3.49) 


3.3  Numerical  Implementation  of  the  Structural  Solver 


3.3.1  Overview 


The  assumed-mode  ordinary  differential  equation  of  motion  ( 3.40)  is  evaluated  numerically  in  FOR¬ 
TRAN  90  using  explicit  integration.  Essentially,  this  process  consists  of  calculating  the  mass, 
damping,  and  stiffness  matrices  and  the  force  vectors  at  each  time  step. 


3.3.2  Matrices 

The  time-invariant  linear  matrices  are  calculated  only  once  using  Simpson’s  ^  rule.  The  forcing 
vectors  and  nonlinear  matrices  are  dependent  upon  time  and,  consequently,  calculated  for  every 
time  step. 

From  a  computational  perspective,  the  triple  integrals  in  the  nonlinear  mass  and  damping 
matrices  are  particularly  expensive.  However,  since  the  deformations  are  expressed  in  terms  of 
separable  variables,  the  time-dependent  contribution  can  be  extracted,  enabling  the  spatial  integrals 
to  be  computed  once  and  accordingly  multiplied  by  the  corresponding  time  coefficients  and  summed 
subsequently. 

For  example,  in  the  triple  integral 


rL  rs  rs 

j  W-(s)w'(s,t)  j  m(s)  J  v'(s,t)Wj(s)dsdsds, 


,  a=l 


.6=1 


The  time  coefficients  and  summations  are  rearranged  to  yield 

^  m  C  rL  rL  rs  \ 

wa(t)vb(t)  W'(s)W^{s)  m(s )  Wl{s)W'(s)dsdsds\ . 
0=11=1  1  Js  J 


(3.50) 


the  deformations  are  expanded  in  terms  of  mode  shapes  and  time  coefficients: 

rL  (  1  )  fL  fS  (  m  1 

l  1  !>.(«»)  f  J  m(s)j  lY.  Vbi^WKs)  >  Wj(s)dsdsds.  (3.51) 


(3.52) 


Furthermore, 

W{,(s)W'(s)ds  (3.53) 

has  a  closed  form  solution,  essentially  reducing  the  calculation  of  the  triple  integral  at  every  time 
step  to  a  once-computed  double  integral  that  undergoes  double  summation  at  subsequent  time  steps. 
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The  remaining  integrals  in  the  nonlinear  matrices  are  calculated  in  a  similar  manner,  reducing  the 
integrals  to  double  summations. 


3.3.3  Forcing  Vectors 

The  forcing  vectors  are  computed  at  every  time  step  using  Simpson’s  ^  rule.  For  the  case  of  p  point 
loads  F;(f),  expressed  in  the  inertial  frame, 


Qu  |  p 

q  =  { Qv\  =  ]rj(s-s*)Ff, 

Qi , 


i= 1 


where  <5(s)  is  the  Dirac  delta  function  ([Habcrman,  2004).  This  is  done  similarly  for 
moment  about  the  £-axis,  leading  to 

fL  p 

{FL}i  =  /  WiQwds  =  YJWi{sj)Fz.{t) 

Jo  j=i 

rL  p 

{FL}l+i  =  WiQvds  =YJWi{sj)FVj(t ) 

J°  3=1 

rL  q 

{FL}l+m+i  =  /  <5>iQ<pds  = 

Jo  ,  , 


(3.54) 
with  the 


(3.55) 


k= 1 


and 


{ FNL)i 

7L 

r  « 

1 

II 

"T3 

3 

o* 

£ 

1 

II 

fsj 

/  Wi(s] 

)w\s ,  t)ds 

JO  Js  j=1 

Jo 

fL  fL 

fS3 

{FNLh+i 

=  -/  WiV'  Quds  =  -Y^FXj(t) 

/  Wi(s] 

)v'(s,  t)ds, 

Jo  Js  -=1 

Jo 

p  i 

{Fnl\i  =  -EE  FXj{t)wa(t) 

3=1  a=1 

p  m 

{FNL}i+i  =  -EE  FXj{t)va(t) 

3= 1  0=1 

where  the  integral  has  a  closed- form  solution. 


f  Wi{s)W'a{s)ds 
Jo 

P  W^W'^ds, 

Jo 


(3.56) 


(3.57) 


3.3.4  Solution 

Letting  the  time  coefficients  be  represented  by 

X  =  {wi,...,wt,  Vi  01,  -  -  - ,  <j>n},  (3.58) 
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equation  (|3.40|)  is  rewritten  such  that 


d_ 

dt 


X 


[Ml  +  Mnl] -1  {Fl  +  Fnl  -  [Kl  +  Knl]  X  -  [CNL]  X} 


or  in  shorthand  as 


(3.59) 


Y  =  F(Y,t). 


(3.60) 


At  each  point  in  time,  (3.59)  is  integrated  using  a  fourth-order  Runge-Kutta  method  (Gerald  and| 
Wheatley|,  |2004|): 


Yn+i  —  Yn  +  -A t  (ki  +  2k2  +  2k3  +  k4) , 
6 


(3.61) 


where 


tn+ 1 

ki 

k2 

k3 

k4 


tn  +  At 

(3.62) 

F  (Yn,  tn) 

(3.63) 

F  ^Yn  +  -Atki ,tn  +  —  At^j 

(3.64) 

F  ^Yn  +  — Atk2,tn  +  —A t'j 

(3.65) 

F  (Yn  +  Atk3,  tn  +  At) . 

(3.66) 

This  provides  the  time  coefficients  and  the  first  derivatives  of  the  time  coefficients  for  each  time 
step.  Consequently,  the  deformation  of  the  elastic  axis  can  be  calculated  with  respect  to  time 
by  summing  the  product  of  the  time  coefficients  and  the  mode  shapes  at  the  position  of  interest. 
Knowing  the  deformed  position  of  each  point  along  the  elastic  axis  as  well  as  < f>,  the  angle  about  £, 
which  the  cross  section  coordinate  system  is  rotated,  enables  the  displacement  of  any  point  on  the 
structure  to  be  calculated. 


3.4  Structural  Solver  Validation 

3.4.1  Method 

The  structural  solver  is  validated  by  comparing  results  with  those  generated  by  Abaqus,  a  commer¬ 
cial  finite  element  software  package.  Validation  consists  of  comparing  frequencies  and  deflections 
obtained  with  the  beam  model  developed  herein  and  the  commercial  finite  element  code  Abaqus. 

3.4.2  Natural  Frequencies 

Three  beams  are  used  as  test  cases  for  comparing  natural  frequencies  obtained  by  using  Abaqus 
and  the  present  beam  model.  Each  of  the  beams  is  modeled  in  Abaqus  using  quadratic  three- 
dimensional  continuum  elements.  Furthermore,  for  each  case,  a  grid  convergence  test  is  done  by 
refining  the  mesh  twice.  The  five  frequencies  corresponding  to  the  first  five  vibrational  modes  are 
compared  with  those  obtained  by  the  linear  contribution  to  the  beam  model.  For  these  comparisons, 
the  beam  model  uses  five  mode  shapes  for  each  of  the  three  independent  degrees  of  freedom. 
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The  linear  contribution  to  the  beam  model  provides  a  system  of  15  coupled  linear  second-order 
differential  equations  for  which  a  sinusoidal  analytical  solution  is  obtainable  for  the  mode  shape 
time  coefficients.  This  solution  provides  the  natural  frequencies  of  the  beam. 


3.4. 2.1  Case  1:  A  Tapered  Beam 


The  first  beam  used  for  frequency  comparison  is  a  tapered,  homogeneous  beam  composed  of  alu¬ 
minum  alloy  6061-T6.  This  beam  choice  provides  an  initial  example  of  a  non-uniform  beam. 
6061-T6  has  a  density  of  2710  kg/m3,  a  Young’s  modulus  of  70  GPa,  and  a  shear  modulus  of  26 
GPa  (Beer  et  al.,  2006|).  The  beam  is  shown  in  Figure  3.4.  It  is  20  cm  long  and  consists  of  a  2  cm  x 
1  cm  cross  section  at  the  fixed  end  and  a  1  cm  x  1  cm  cross  section  at  the  free  end.  The  properties 
of  the  tapered  beam  are  listed  in  Table  3.1.  For  a  rectangular  cross  section  with  dimensions  a  x 


b,  where  a  is  the  greater  of  the  two,  the  torsion  constant  can  be  approximated  as  follows  (Roark 


et  al.|,  |2002|): 


K  =  abA 


1  b  (  bA 

-  -  0.21-  1 - 7 

3  a  V  12a4 


(3.67) 


Figure  3.4:  Tapered  beam  with  100x5x6  mesh. 


Table  3.2  shows  the  mesh  sizes  for  the  three  finite  element  meshes.  For  each  of  the  three 


meshes,  the  first  five  frequencies  are  tabulated  in  Table  3.3.  The  error  is  measured  by  taking  the 
absolute  value  of  the  difference  between  the  frequencies  of  the  finest  finite  element  mesh  and  the 
assumed-mode  solution  and  dividing  by  that  of  the  finite  element  mesh. 

The  assumed-mode  solution  favorably  agrees  with  the  finite  element  solution,  providing  little 
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Table  3.1:  Properties  of  tapered  beam. 


Parameter 
m(s ) 

jds) 

Ms) 

Jds) 

jr/dS) 

Dds ) 


Dv(s) 
Dds) 
Dvds ) 


Value 


0.542  -  1.355s  kg/m 

(0.22583  -  1.4679s  +  3.3875s2  -  2.8229s3)  x  10"4 
(0.18067  -  1.3550s  +  3.3875s2  -  2.8229s3)  x  10"4 
(0.045167-  0.11292s)  x  10"4  kg  m2/m 


0  kg  nr2 /nr 
173.33  + 


9.1xlO-10 


1.092 


(0.02— 0.05s)5  (0.02— 0.05s) 

loo  00n  (2.275x10_9s)  ,  (2.73s) 

400.  -MS  (0.02— 0.05s)5 


N  m2 


(0.02— 0.05s) 

466.67  -  3500s  +  8750s2  -  7291.7s3  N  m2 
116.67-  291.67s  N  m2 


kg  m2/m 
kg  m2/m 


0  N  m2 


Table  3.2:  Finite  element  mesh  sizes  of  tapered  beam. 


Mesh  Number 

Elements  Along  x,  y,  z 

Total  Elements 

1 

50  x  3  x  3 

450 

2 

67  x  3  x  5 

1005 

3 

100  x  5  x  6 

3000 

Table  3.3:  First  five  frequencies  of  tapered  beam. 


Mode 

Mesh  1  [Hz] 

Abaqus 
Mesh  2  [Hz] 

Mesh  3  [Hz] 

Modal 

Representation  [Hz] 

Error 

1 

252.75 

252.71 

252.64 

251.89 

0.297% 

2 

445.46 

445.38 

445.30 

446.20 

0.202% 

3 

1362.1 

1361.8 

1361.4 

1370.8 

0.686% 

4 

2081.4 

2080.8 

2080.4 

2130.5 

2.41% 

5 

3594.1 

3593.4 

3592.1 

3671.3 

2.16% 
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error. 


3. 4. 2. 2  Case  2:  A  Twisted  Beam 


A  twisted,  homogeneous  beam  composed  of  aluminum  alloy  6061-T6  is  the  second  beam  used  for 
frequency  comparison.  This  beam  is  chosen  as  it  provides  an  example  of  a  non-uniform  beam  for 
which  the  principal  mass  and  bending  axes  vary  with  respect  to  the  reference  coordinate  system. 
The  beam  is  shown  in  Figure  |3.5|.  It  is  10  m  long  and  consists  of  a  1  m  x  0.5  m  cross  section  that  is 
unrotated  at  the  fixed  end  and  linear  rotated  such  that  it  is  at  a  45  angle  about  the  centroidal  axis 
at  the  free  end.  The  properties  of  the  twisted  beam  are  presented  in  Table  [3.4|.  For  a  rectangular 
cross  section  with  an  aspect  ratio  of  2,  the  torsion  constant  K  is  0.229a63,  where  a  and  b  are 
respectively  the  lengths  of  the  longer  and  shorter  sides  (Ugural  and  Fenster,  1981|). 


Figure  3.5:  Twisted  beam  with  100x5x10  mesh. 


As  with  the  tapered  beam,  Table  3.5  provides  the  mesh  sizes  for  the  three  finite  element  meshes 
For  each  of  the  three  meshes,  the  first  five  frequencies  and  the  error  are  tabulated  in  Table 
The  assumed-mode  solution  agrees  nicely  with  the  finite  element  solution. 


3. 4. 2. 3  Case  3:  A  Composite  Beam 

The  final  beam  used  for  frequency  analysis  is  a  uniform  composite  beam  consisting  of  aluminum 
alloy  6061-T6  and  steel.  Steel  has  a  density  of  7860  kg/m3,  a  Young’s  modulus  of  200  GPa,  and  a 
modulus  of  rigidity  of  77.2  GPa  (|Beer  et  ah ,  [2006).  Three  quarters  of  the  beam  are  aluminum,  while 
the  remaining  quarter  is  steel.  The  beam  is  shown  in  Figure  3.6.  This  beam  is  selected  because 
it  provides  an  example  of  a  beam  for  which  the  mass  and  elastic  centers  of  each  cross  section  are 
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Table  3.4:  Properties  of  Twisted  Beam 


Parameter 

Value 

m(s) 

1355  kg/m 

k(s) 

141.15  kg  m2/m 

Ms) 

70.573  +  42.344  cos  kg  m2/m 

k(s) 

70.573  -  42.344  cos  g  kg  m2/m 

3vds ) 

—42.344  sin  kg  m2/m 

D&) 

7.4425  x  108  N  m2 

Dds ) 

(1.8229  +  1.0938  cos  §f )  x  109  N  m2 

D((s) 

(1.8229  -  1.0938  cos  x  109  N  m2 

Dv  c(a) 

-1.0938  cos  ^  x  109  N  m2 

Table  3.5:  Finite  element  mesh  sizes  of  twisted  beam. 


Mesh  Number 

Elements  Along  x,  y,  z 

Total  Elements 

1 

50  x  2  x  5 

500 

2 

67  x  3  x  7 

1407 

3 

100  x  5  x  10 

5000 

Table  3.6:  First  five  frequencies  of  twisted  beam. 


Mode 

Mesh  1  [Hz] 

Abaqus 
Mesh  2  [Hz] 

Mesh  3  [Hz] 

Modal 

Representation  [Hz] 

Error 

1 

4.1440 

4.1423 

4.1410 

4.1291 

0.288% 

2 

8.0006 

7.9981 

7.9961 

8.0229 

0.335% 

3 

26.265 

26.253 

26.245 

26.483 

0.907% 

4 

46.756 

46.737 

46.723 

48.425 

3.64% 

5 

58.102 

58.041 

58.024 

57.407 

1.06% 

49 


offset  in  both  dimensions.  Furthermore,  the  principal  bending  and  mass  axes  are  neither  parallel 
mutually  nor  to  the  reference  coordinate  system.  The  beam  is  1  nr  long  with  a  5  cm  X  3  cm  cross 
section.  The  properties  of  the  composite  beam  are  presented  in  Table  |3.7|.  The  torsional  stiffness 
is  approximated  by  multiplying  an  area-weighted  average  of  the  modulus  of  rigidity  by  the  result 
of  (|3.67)  for  the  cross  section. 


Figure  3.6:  Composite  beam. 


Table  |3.q  tabulates  the  mesh  sizes  for  the  three  finite  element  meshes,  and  Table  T9  lists  the 
first  five  frequencies  and  the  error.  The  frequencies  calculated  by  the  beam  modal  representation 
provided  accurate  results  while  reducing  the  computational  time  by  approximately  two  orders  of 
magnitude. 
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Table  3.7:  Properties  of  composite  beam. 


Parameter 

Value 

ev(s) 

3.7524  x  10-5  m 

ec{s) 

-6.2539  x  10“5  m 

m(s) 

5.9963  kg/rn 

jds) 

1.5668  x  10-3  kg  m2/m 

1.1521  x  10-3  kg  m2/m 

Jds) 

0.41474  x  10-3  kg  m2/m 

jr/dS) 

0.12276  x  10-3  kg  m2/m 

D^s) 

1.0931  x  104  N  m2 

Dds) 

2.9616  x  104  N  m2 

DC(s) 

1.0662  x  104  N  m2 

Dvds ) 

0.31212  x  104  N  m2 

Table  3.8:  Finite  element  mesh  sizes  of  composite  beam. 


Mesh  Number 

Elements  Along  x,  y,  z 

Total  Elements 

1 

77  x  2  x  4 

616 

2 

100  x  4  x  6 

2400 

3 

143  x  4  x  8 

4576 

Table  3.9:  First  five  frequencies  of  composite  beam. 


Mode 

Mesh  1  [Hz] 

Abaqus 
Mesh  2  [Hz] 

Mesh  3  [Hz] 

Modal 

Representation  [Hz] 

Error 

1 

23.099 

23.091 

23.089 

23.036 

0.230% 

2 

39.675 

39.668 

39.665 

39.662 

0.00756% 

3 

144.18 

144.13 

144.12 

144.30 

0.125% 

4 

245.73 

245.68 

245.66 

248.21 

1.04% 

5 

401.19 

401.01 

400.98 

403.61 

0.656% 
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3.5  Finite  Element  Method  Code  for  Plates 


Although  the  nonlinear  beam  element  presented  above  excels  on  both  accuracy  and  efficiency,  not 
all  structures  can  be  modeled  by  beams.  For  this  reason,  and  because  we  prefer  to  have  access  to 
the  source  code  rather  than  using  canned  commercial  codes,  during  this  project  we  also  developed 
a  finite  element  method  for  plates. 

The  Bogner-Fox-Schmit  (BFS)  rectangular  conforming  element  was  considered  for  rectangular 
and  triangular  finite  elements  used  herein  (Bogner  et  alj,  1966|).  The  FEM  employed  three- node 
triangular  Mindlin  (MIN3)  plate  element  with  an  improved  transverse  shear  term  ([Tessler  and 
Hughes,  1985).  The  nodal  degree  of  freedom  (DOF)  of  each  BFS  element  was  composed  of  16 
bending  and  8  in-plane  DOF.  The  MIN3  element  has  a  total  of  15  DOF,  5  at  each  apex  node.  The 
bending  node  DOF  comprises  transverse  displacements  and  normal  rotations,  and  the  in-plane  node 
DOF  comprises  in-plane  displacements.  The  elements  have  been  studied  for  similar  cases,  including 
composite  panels  under  the  combined  aerodynamic  and  acoustic  pressures  based  on  modal  analysis. 
Further,  the  MIN3  triangular  element  was  used  to  study  LCO  of  arbitrary  laminated  anisotropic 
composite  rectangular  plates.  The  FEM  used  the  BFS  and  Min3  elements.  The  Generic  Transport 
Wing  (GTW),  as  described  by  Edwards  et  al.  (2009),  was  examined  in  detail. 

The  Generic  Transport  Wing  (GTW)  provided  us  a  platform  for  computational  fluid  dynam¬ 
ics/computational  structural  mechanics  experimental  validation.  The  GTW  was  used  in  legacy 
experiments  within  the  NASA  LaRC  TDT  Facility.  Features  of  the  GTW  include:  a  structure 
of  solid  aluminum  with  balsa  skin,  LCO  responses  exhibited  within  the  transonic  regime,  LCOs 
exhibited  with  a  winglet  present  and  noticeable  differences  in  vibration  characteristics  from  the 
baseline,  and  repeatable  measurements  of  LCO.  The  GTW  geometry  is  shown  in  the  Fig.  |3.7|. 

Indeed,  the  LCO  was  present  on  the  GTW  with  winglet  (as  well  as  the  other  configurations), 
and  the  winglet  made  noticeable  differences  in  the  vibration  characteristics.  The  plate-based  FEM 
model  we  developed  was  used  to  replicate  the  results  published  by  Edwards  et  al.  ( 2009|) .  The 


vibration  characteristics  are  given  in  Table  |3.10j  and  the  first  three  modes  are  illustrated  in  Fig. 


Table  3.10:  First  four  frequencies  of  generic  transport  wing. 


Mode 

Measured 

[Hz] 

FEM 

[Hz] 

Error 

[%] 

1st  bending 

4.08 

4.07 

0.245 

2nd  bending 

14.0 

14.5 

3.571 

1st  torsion 

31.5 

29.3 

6.984 

2nd  torsion 

58.1 

57.4 

1.205 
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Generic  Transport  Wing  (GTW) 


Figure  3.7:  Generic  transport  wing. 
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l“  Bending  Mode 


2"d  Bending  Mode 


1  *'  Torsion  Mode 


Figure  3.8: 


Generic  transport  wing:  first  three  modes. 
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Chapter  4 

Continuation  Methods 


Since  the  design  space  of  new  air  vehicles  is  vast,  determination  of  nonlinear  responses,  vehicle  trim 
states  and  vehicle  stability  must  be  addressed  in  an  efficient,  methodical  manner.  Numerical  Con¬ 
tinuation  will  reveal  various  bifurcation  characteristics  including  both  subcritical  and  supercritical 
bifurcations.  The  former  case  is  of  significant  interest  because  the  subcritical  bifurcations  indicate 
limit  cycle  oscillations  below  the  linear  stability  boundary.  Further,  subcritical  bifurcations  depend 
upon  disturbances,  and  have  a  hysteresis  between  onset  and  recovery  speeds.  This  behavior  cannot 
be  studied  with  linear  methods.  Subcritical  bifurcations  are  anticipated  for  highly  deformable,  high 
aspect  ratio  wings  and  may  reveal  adverse  vehicle  performance. 

Numerical  continuation  methods  offer  a  robust  and  efficient  solution  to  these  needs.  In  our 
opinion,  one  limitation  is  the  ability  to  properly  represent  the  aerodynamic  loads  as  appropriate 
for  the  regime.  However,  with  the  inclusion  of  the  variable-fidelity  reduced-order  model,  the  Con¬ 
tinuation  Method  permits  direct  analysis  of  the  system  without  requiring  explicit  linearization  or 
compromise  on  a  physical  description  of  the  flow.  Thus,  the  trim  states  of  the  nonlinear  system  and 
the  stability  characteristics  of  the  air  vehicle  can  be  evaluated  in  a  continuous  manner  as  system 
parameters  are  varied.  Continuation  Methods  are  direct  methods  for  examining  the  bifurcation 
characteristics  of  the  nonlinear  equations. 

We  use  the  software  tool  AUTO  for  continuation  and  bifurcation  analyzes  of  the  nonlinear 
ordinary  differential  equations  flDoedel,  1981|).  AUTO  automates  the  computation  of  solutions  of 
this  parameter-dependent  system  of  equations.  The  system  possesses  multiple  solutions,  and  it 
is  valuable  to  compute  the  set  of  solutions  and  search  for  those  solutions  with  specific  desirable 
properties  as  a  system  parameter  is  varied.  This  solution  set  forms  a  bifurcation  diagram,  i.e. , 
a  smooth  curve  (or  surface)  representing  the  solution  for  the  varying  system  parameter.  Tools 
such  as  AUTO  facilitate  parametric  studies  and  minimize  time  consuming  numerical  simulations. 
The  computation  of  such  bifurcation  diagrams  and  associated  singularities  is  captured  within  the 
domain  of  the  numerical  continuation  algorithm.  AUTO  has  the  capability  to  compute  periodic 
solutions  such  as  the  limit  cycle  oscillations  (LCOs).  In  addition,  AUTO  determines  the  stability 
of  the  steady  state  solutions.  AUTO  permits  a  two-parameter  continuation  solution  of  Hopf  and 
other  bifurcation  points.  A  detailed  description  on  the  features  of  AUTO  can  be  found  in  (pocdel 
20021). 


et  al 


The  method  permits  direct  eigen-analysis  of  the  system  without  resorting  to  explicit  lineariza¬ 
tion.  Thus,  the  nonlinear  trim  states  and  their  stability  characteristic  can  be  evaluated  in  a  contin¬ 
uous  manner  as  a  single  system  parameter  is  varied  using  the  exact  nonlinear  expressions.  Eigen¬ 
values  obtained  by  the  method  are  used  to  ascertain  the  extent  of  coupling  between  rigid  body 
and  flexible  body  modes.  Various  bifurcation  points  of  the  system  are  found  and  permit  evaluation 
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of  stability  boundaries.  Also,  the  method  is  used  to  examine  the  presence  of  LCOs.  In  such  a 
situation,  the  dynamic  characteristic  of  the  system  is  dependent  on  initial  condition. 

We  have  applied  the  continuation  method  to  the  nonlinear  aeroelastic  responses  for  vehicles 
with  moderate-to-high  deformations.  The  wing  was  modeled  using  nonlinear  structural  equations 
of  motion  for  a  cantilevered  beam  with  in-plane,  out-of-plane,  and  torsional  couplings.  In  the  form 
of  a  bifurcation  diagram,  Fig.  4.1  is  an  example  of  a  bifurcation  diagram  created  as  part  of  the 
continuation  tool.  The  tool  permits  nonlinear  analysis  of  the  set  of  governing  equations.  The 
approach  detects  the  bifurcation  points,  determines  equilibria  and  stability  states,  determines  am¬ 
plitude  of  periodic  solutions  (e.g.,  stable  LCOs),  and  rapidly  identifies  character  of  design.  We  have 
shown  that  it  duplicates  time-domain  simulation,  but  such  exhaustive  simulation  to  “bracket”  the 
stability  boundaries  is  avoided. 


Velocity  V  [ft/s] 


Figure  4.1:  Example  bifurcation  diagram  showing  the  effect  of  trim  constraint. 


Figure  LI  shows  amplitudes  of  the  out-of-plane  modal  coordinate,  equilibria,  and  limit  cycles 
as  the  freestream  velocity  is  varied.  Bifurcation  occurs  at  the  critical  velocity.  From  the  bifurcation 
point,  the  thin  solid  line  represents  a  branch  of  stable  equilibria,  the  thin  dashed  line  represents 
a  branch  of  unstable  equilibria,  and  the  thick  solid  line  represents  a  branch  of  stable  LCOs.  An 
important  observation  is  that  the  trim  model  shows  a  change  in  amplitude  in  both  equilibria  and 
LCOs.  In  addition,  the  continuation  tool  allows  one  to  explore  the  parameter  space  without  missing 
details. 

Figure  |4.2|  shows  results  for  the  wing  with  a  trim  constraint  such  that  the  wing  possesses  a  root 
angle.  Under  flight  conditions  requiring  a  root  angle  of  attack  to  satisfy  the  trim  constraint,  the 
static  deformation  and  stability  characteristics  change.  The  trim  condition  leads  to  a  statically 
deformed  shape  of  the  wing  that  must  vary  with  the  freestream  velocity  and,  as  one  might  expect, 
static  deformation  of  the  wing  increases  with  increasing  freestream  velocity.  The  system  loses 
stability  at  a  velocity  of  approximately  124  ft/sec. 

The  left  view  of  Fig.  |4.2|  shows  a  bifurcation  diagram  with  the  supercritical  branch  which 
describes  the  response  at  velocities  higher  than  the  bifurcation  point.  The  nature  of  bifurcation 
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V  •  I  ««  If  f  V  I  If 


AR  =  33  AR  =  40 

Figure  4.2:  Sensitivity  to  aspect  ratio:  transition  to  subcritical  from  supercritical  behavior  at  angle 
of  attach  2.8  deg  and  in-plane  to  out-of-plane  stiffness  ratio,  (3n= 4. 


depends  on  the  angle  of  attack  imposed  by  the  trim  condition  as  well  as  system  design  properties 
such  as  the  in-plane  to  out-of-plane  stiffness  ratios  (/3n).  hi  addition,  the  right  view  indicates 
the  vehicle  transitions  to  a  subcritical  branch  at  a  higher  aspect  ratios.  In  linear  analysis,  the 
bifurcation  point  is  the  flutter  velocity,  and  is  the  only  response  predicted  by  such  analysis.  The 
supercritical  branch,  as  well  as  the  stable  branch  (solid  line)  in  the  right  view,  describes  limit  cycle 
behavior.  The  nonlinear  dynamics  of  this  system  are  governed  by  the  existence  of  a  disturbance. 
With  a  disturbance,  the  system  becomes  entrained  in  the  LCO  with  amplitude  as  shown  by  the 
thick  solid  line. 

In  summary,  the  Continuation  Method  provides  an  efficient  method  of  nonlinear  analysis  for 
the  aeroelastic  system.  The  method  eliminates  “bracketing”  approaches  inherent  of  time  domain 
simulation  of  nonlinear  equations.  Continuation  detects  bifurcations,  defines  equilibria  and  stability, 
determines  amplitude  of  LCOs,  characteristics,  and  identifies  design  characteristics  in  a  rapid, 
comprehensive  manner. 
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Chapter  5 

A  Novel  Flow  Solver  Grid  Generator 
for  Aircraft  with  Extreme  Wing 
Deformation 


Grid  generation  is  an  important  part  of  computational  fluid  dynamics,  which  greatly  affects  the 
accuracy  of  the  numerical  results.  Grid  quality  has  an  even  stronger  effect  on  solution  accuracy  when 
simulating  unsteady  flows,  and  especially  when  the  mesh  deforms,  like  in  aeroelasticity  simulations. 
This  chapter  presents  a  new  grid  generation  algorithm  for  aircraft  with  extreme  wing  deformation. 
The  next  section  briefly  presents  background  information  on  the  grid  generator  that  was  enhanced 
during  this  project.  The  following  section  presents  the  conformal  mappings  proposed  to  enhance 
the  existing  grid  generator.  The  last  section  compares  the  grid  qualities  of  the  two  grid  generators. 


5.1  Background  information 


The  new  grid  generation  algorithm  proposed  herein  was  implemented  into  an  existing  mesh  gen¬ 
eration  code  flCizmas  and  Gargoloff,  200S|).  This  existing  grid  generation  code  discretized  the 
computational  domain  of  the  flow  solver  using  a  hybrid  grid  consisting  of  hexahedra  and  triangular 
prisms.  Because  of  the  wing-fuselage  configuration,  it  was  possible  to  divide  the  computational 
domain  into  layers  that  were  topologically  identical,  as  shown  in  Fig.  5.1.  Therefore,  a  2.5-D  grid 


could  be  used  to  discretize  the  3-D  domain,  which  simplified  the  grid  generation  algorithm. 

The  topologically  identical  nodes  of  adjacent  layers  were  interconnected  to  generate  the  volume 
elements  that  were  either  triangular  prisms  or  hexahedra.  To  better  control  mesh  size  near  the  wing, 
each  layer  had  a  structured  O-grid  around  the  wing,  and  an  unstructured  grid  at  the  exterior  of 
the  O-grid.  The  O-grid  allowed  clustering  cells  in  the  direction  normal  to  the  structure  to  properly 
capture  boundary  layer  effects.  The  unstructured  grid  was  flexible  in  filling  the  rest  of  the  domain. 

The  O-grid  algorithm  previously  used  (|Cizmas  and  Gargolof!|,  2008|)  was  based  on  a  Poisson 
solver  which  did  not  enforce  grid  orthogonality.  As  a  result,  the  quality  of  the  O-grid  deteriorated 
near  the  leading  and  trailing  edges  of  the  wing.  To  alleviate  this  problem,  a  conformal  mapping 
algorithm  was  developed  herein,  which  assured  grid  orthogonality  over  the  entire  O-grid.  Grid 
orthogonality  is  also  beneficial  for  the  cases  when  wing  chord- wise  bending  is  significant.  In  these 
cases,  a  good  quality  grid  based  on  the  Poisson  solver  is  difficult  to  obtain  for  all  chord- wise  bending 
configurations. 
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(a) 


(b) 


Figure  5.1:  Grid  generator:  (a)  detail  of  topologically  identical  layers  along  the  wing,  and  (b) 
aircraft  configuration. 


5.2  Conformal  Mapping 


The  conformal  mapping  used  herein  to  generate  the  O-grid  is  an  extension  of  the  Carafoli  method 
for  determining  the  velocity  and  pressure  distribution  on  airfoil  cascades,  assuming  the  flow  is 
inviscid,  incompressible  and  potential  (Berbente  and  Cizmas,  1990|).  This  conformal  mapping  is 
used  for  two  purposes:  (1)  to  generate  the  grid  around  the  wing,  and  (2)  to  provide  a  start-up  value 
for  the  velocity  and  pressure  field  in  the  O-grid. 

The  conformal  mapping  proposed  herein  consisted  of  three  transformations  that  mapped  the 
exterior  of  a  unit  circle  onto  either  an  isolated  airfoil  or  an  airfoil  in  a  cascade.  It  is  known  that 
when  the  distance  between  airfoils  is  larger  than  approximately  three  chords,  the  effect  of  the 
neighboring  airfoil  is  negligible.  As  a  result,  the  user  has  the  option  to  generate  a  grid  for  either  an 
isolated  airfoil  or  a  cascade  airfoil  by  selecting  the  distance  between  airfoils.  Therefore,  this  grid 
generation  algorithm  can  also  be  used  for  turbomachinery  flows. 


The  first  conformal  mapping  of  the  sequence  is  the  Weinig-Manea  transformation  (Berbente 
and  Cizmas,  1990|) 


z  = 


2vr 


=  —  I  eiX  In 


(5.1) 


where  t  was  the  cascade  pitch,  A  was  the  stagger  angle,  (  was  the  complex  variable  of  the  mapped 
circle  domain,  2  was  the  complex  variable  of  the  airfoil  domain,  and  R  and  B  were  constants  that 


were  determined  as  a  function  of  the  solidity  and  stagger  angle.  Function  (5.1)  maps  the  exterior 
of  the  unit  circle  into  a  cascade  of  flat  plates.  Using  (|5.1|)  with  the  translation  and  homothety  of 
the  unit  circle  generates  cascade  airfoils  with  various  thicknesses  and  curvatures. 

To  generate  the  grid  and  find  the  velocity  field  for  a  given  airfoil,  a  conformal  mapping  similar 


to  (5.1)  would  be  needed,  which  maps  the  unit  circle  onto  the  wing  airfoil.  To  achieve  this  goal, 
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a  series  of  conformal  mappings  are  determined  by  going  backward,  that  is,  from  the  airfoil  to  the 
unit  circle. 


<S,  z2  tf.z, 


The  sequence  of  mappings,  shown  in  Fig.  5.2,  consists  of:  (1)  conformal  mapping  2(22)  given 
by  (p|)  where  C  was  substituted  by  Z2,  (2)  conformal  mapping 


22 


eiT  +  z0 


(5.2) 


which  represents  the  mapping  of  a  circle  to  an  ellipse,  where  r  is  the  angle  of  the  axes  of  the  ellipse, 
c2  =  ( b 2  —  a2)/ 4  with  a  and  b  being  the  axes  of  the  ellipse,  and  (3)  conformal  mapping 


m— 2 

21  =  ]T  C.jCj,  c.j  €  C,  (5.3) 

3=~ 1 

where  the  constants  C-j  were  determined  given  a  contour  C  in  plane  cto. 

The  implementation  of  the  numerical  algorithm  that  calculates  the  conformal  mapping  is  briefly 
described  herein.  For  a  given  airfoil  shape,  the  Newton-Raphson  method  was  used  to  determine 
the  shape  of  the  contour  E  in  plane  $ .  The  closest  ellipse  to  contour  E  was  then  obtained  by 
determining  zo,  a,  b  and  r  using  the  least  squares  method.  The  coordinates  of  contour  C  in  plane 
cto  were  determined  using  a  Newton-Raphson  method  applied  to  the  conformal  mapping  Q5.2|).  The 
complex  coefficients  CTj  were  determined  by  rewriting  fl5.3|)  in  polar  coordinates 

m— 2 

z i  =  ^2  (A -j  +  iB-j)r~^ (cos  j<p  —  isinjtp), 
j=- 1 


dividing  the  unit  circle  K  in  2  m  equal  segments  (m  =  2k.  k  £  N)  and  using  the  orthogonality  of 
trigonometric  functions  with  discrete  arguments. 

Once  the  conformal  mappings  22(21)  and  21(C)  are  known,  the  conformal  mapping  from  the 
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unit  circle  to  the  airfoil  was  calculated  as 


zp  =  z  o  z2  °  z\ 

Figure  [O]  shows  a  detail  of  the  structured  and  unstructured  grids  as  seen  along  the  y-axis.  The 
structured  O-grid  was  generated  using  the  conformal  mapping  approach  presented  above. 


Figure  5.3:  Conformal  mapping  O-grid  surrounded  by  unstructured  grid. 


5.3  Grid  Quality 


The  quality  of  the  unstructured  mesh  was  investigated  in  (|Cizmas  and  Gargolofl ,  2008|) .  To  evaluate 
the  quality  of  the  structured  mesh,  an  indicator  based  on  angles  is  proposed  herein.  This  indicator 
assessed  the  quality  of  the  layers  perpendicular  to  the  wing  span  direction.  Therefore,  the  three- 
dimensional  grid  quality  was  reduced  to  a  two-dimensional  grid  quality  problem. 

Based  on  angles,  the  cell  with  the  highest  quality  is  a  rectangle.  Consequently,  the  quality 
indicator  was  defined  as  q  =  0.25  ^T=1  4  sinai,  where  is  the  angle  at  corner  i.  As  a  result,  the 
mesh  quality  improved  as  the  indicator  got  closer  to  1,  the  maximum  value  of  q. 

Figure  |5.4|  shows  the  contour  plots  of  the  mesh  quality  indicator  of  an  O-grid  generated  using  the 
Poisson  solver  and  the  conformal  mapping.  The  grid  generated  using  the  Poisson  solver  had  6250 
nodes  while  the  grid  generated  using  the  conformal  mapping  had  6400  nodes.  The  average  angle 
for  the  Poisson  solver  was  63  deg,  while  for  the  conformal  mapping  was  81  deg .  Note  that  although 
the  conformal  mapping  generates  an  orthogonal  grid,  the  edges  of  the  cells  are  approximated  by 


straight  lines  and  therefore  the  average  angle  was  not  90  deg.  Figure  5.4  shows  that  the  grid  quality 
of  the  conformal  mapping  grid  is  superior  to  that  of  the  Poisson  solver. 

To  test  the  conformal  mapping  algorithm  for  extreme  cases  of  cross-wise  bending,  the  wing 
deformation  was  amplified  as  shown  in  Fig.  |5.5|.  As  expected,  the  quality  of  the  grid  generated 
using  conformal  mapping  was  good,  in  spite  of  the  large  cross- wise  deformation. 

The  grid  generation  algorithm  was  designed  such  that  the  grid  connectivity  remained  unchanged, 
while  the  mesh  deformed  as  the  wing  bent  and  twisted  under  aerodynamic  loads.  The  algorithm 
for  grid  deformation  was  applied  in  two  steps:  (1)  grid  translation,  and  (2)  grid  rotation  about 
the  three  axes  (Cizrnas  and  Gargolofl,  2008).  The  O-grid  layers  deformed  as  the  wing  cross- 
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Figure  5.4:  Grid  quality  based  on  angle:  (a)  Poisson  solver,  and  (b)  conformal  mapping. 


section  deformed.  The  unstructured  grid  was  deformed  in  the  xz- plane  using  a  spring  analogy 
technique  (Batina,  1990|).  The  nodes  of  each  layer  were  also  translated  in  the  y-direction  according 
to  the  wing  deformation. 

The  grid  was  rotated  and  stretched  about  the  x-,  y-,  and  2-axes.  The  numerical  values  of 
the  rotation  angles  9X,  6y ,  and  9Z  were  calculated  using  data  provided  by  the  structural  solver 
(Cizrnas  and  Gargoloff,  200S|).  Figure  |5T|  shows  the  undeformed  and  deformed  grids  for  the  Goland 
wing  (Eastep  and  Olsen,  |1980|;  |Snyder  et  al.,  2003|) .  In  this  case,  the  deformation  at  the  tip  was 
approximately  60%  of  the  wing  semi-span  length,  illustrating  that  the  grid  deformation  algorithm 
is  robust  and  precludes  the  need  for  regriding,  even  at  large  deformations. 

This  grid  generation  and  deformation  algorithm  allowed  the  reduction  of  the  computational 
cost  of  grid  generation  by  a  factor  of  3.6  or  more  (Cizrnas  and  Gargoloff,  2008|),  while  maintaining 
the  grid  topology  unchanged.  Having  an  unchanged  topology,  while  the  grid  deformed,  simplified 
both  the  parallel  and  multigrid  algorithms. 
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Figure  5.5:  Hybrid  grid  around  the  wing  with  O- 
grid  generated  using  conformal  mapping,  assuming 
exaggerated  chord-wise  deformation. 


Figure  5.6:  Isometric  projection  of  deformed 
and  undeformed  grids  of  the  Goland  wing  -  tip 
deformation  of  60%  of  wing  semi-span  (Piz- 
mas  and  Gargolofl,  2008|). 
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Appendix  A 

Nonlinear  Beam  -  Hamilton’s 
Principle 


The  extended  form  of  Hamilton’s  principle  is 
rt  2 
It! 


51  =  6 


rt2 

/  (T  +  W) 

J  t-\ 


dt  =  0. 


Noting  that 


C  =  T  -  V 
W=  Wc  +  Wnc 
5V  =  -  5WC, 


Hamilton’s  principle  can  be  rewritten  as  Meirovitch]  (|1967) 


51  =  5 


f  (C  +  Wnc)  dt  —  0, 

Jti 


where 


SWnc  =  5Wb  +  [  ( Qu5u  +  Qv6v  +  Qw6w  +  Q^Scj))  ds. 

Jo 


This  leads  to 


rt2 


51  = 


5L  +  5Wb  + 


\L  (Qt 
Jo 


,5u  +  Qv5v  +  Qw5w  + 


ds  >  dt  =  0, 


which,  upon  accounting  for  the  inextensionality  constraint,  becomes 


61  = 


{/  (^51  +  5  -A  ^1  -  (l  +  u')2  +  v'2  +  w'2^j  'j  ds 
6Wb  +  J  (Qu5u  +  Qv6v  +  Qw6w  +  Q(f>5(j))ds^dt  =  0. 


(A.l) 


(A.2) 

(A.3) 

(A.4) 


(A.5) 


(A.6) 


(A.7) 


(A.8) 
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The  variations  of  the  cross-sectional  Lagrangian,  0,  and  p  are 


..  di.  di.  di.  di,  di  .  •  ai 

61=  -wr8u  + —8v  + —Sw  + —Sip  + —^8p  + —8p> 

du  dv  dw  dp >  Op  dp' 


dv 


dw 

dl  dl  ■  dl  ,  dl  dl  ■  dl  , 

dO  do  90'  dp  dip  dp' 

do 


80  = 
8p  = 


,  80  ,  dO  , 

-8u  +  —Sv  +  -w—,Sw 
du  dv  dw 


dip 


8u' 


dip 


8v' . 


du'  dv' 

The  variation  of  the  cross-sectional  Lagrangian  is  then  written  as 

dl  dl  dl  dl  f  dip  /  dip  , 

du  dv  dw  dip  \du'  dv' 

dl  (  d2p  ,  dip  ,  d2p  ,  dip  , 

+  dpi  \dw^5u  +  w8u  +  di^'8v  +  w8v 
dl  (  d2ip  ,  dip  „  d2p  ,  dip  „ 

du'  dsdv  dv 


d^'  \dsdu' 

dl  (dO  !  dO  .  ,  dO  ,  . 

+  m\mSu  +  wSv  +  du 7Sw  |w 


+  -r- 


dl  (  d20 


dO 


d20 


dO  \dtdu'8u  du'88  dtdv’ 


dO 


d20 


5v>  +  w8*'  +  dtdw' 


r  /  8)0  , 

dw  +  — — -dw 
dw' 


dl  (  d20 


..  ,  86  „ 

8u  +  — — -8u  + 


dO'  \dsdu'  U  du' 
dl  dl  •  dl 

+ 

dp  dp  dp' 


d20 


dO 


dsdv'8v  dv'8  dsdw 


d20  .  ,  dO  .  „ 

8w  +  77 — 7  8w 

dw’ 


(A.9) 

(A.10) 

(A. 11) 


(A- 12) 


Expanding  (|A.8|)  and  performing  integration  by  parts  to  express  the  integrand  in  terms  of  the 
variations  8u,  8v,  8w ,  and  Sp  yields 
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51  = 


Letting 


rt2  rl 

Ju  Jo 


+ 


+ 


+ 


+ 


+ 


+ 


+ 


+ 


d3l 


d2l  dip 
dsdip  du' 
dip  d2l 


dl  d2ip  d3l  dip  d2l  d2ip 
dip  dsdu'  +  dtdsdip  du'  +  dtdip  dsdu’ 


5u 


ds2dip '  du' 
d3l  d9 
dtdsdO  du' 
d2l  dip 
dsdip  dv' 
d3l  dip 
ds2dip'  dv' 
d3l  d6 
dtdsdO  dv' 
d2l  dd 
dsdQ  dw' 
d2l  d20 
dtdO  dsdw' 
dl_  _  d2l 
dcp 


+ 


d2ip  d2l  dd  dl  d26 


+ 


dsdip'  dsdu' 
d2l  d26 


+ 


dsdO  du' 
d3l  d6 
ds2dO'  du' 
d3l  dip 


5u 


+ 


dtdO  dsdu' 
dl  d2ip 

dip  dsdu1  +  dtdsdip  dv' 
d2l  d2ip  d2l  dO 


+ 


+ 


dO  dsdu' 
d2l  d20 
dsdO'  dsdu' 
d2l  d2ip 
dtdip  dsdv' 
dl  d20 


5u 


5v 


+ 


dsdip'  dsdv' 
d2l  d20 


dtdO  dsdv' 
dl  d20 


+ 


+ 


+ 


dO  dsdw' 
d3l  dO 
ds2d0'  dw' 
d2l 


dsdO  dv' 
d3l  dO 
ds2d0'  dv' 
d3l  dO 
dtdsdO  dw' 
d2l  d20 


5v 


dO  dsdv' 
d2l  d20 
dsdO'  dsdv' 

5w 


+ 


dv 


+ 


dtdcp  dsdcp' 


5<p  — 


dsdO'  dsdw' 
d2l  d2l 


-5u  — 


5w 


5v  — 


d2l 


5w 


dtdu  dtdv  dtdw 
+  [A7  (l  +  u')  +  A u"]  5u  +  [AV  +  Xv"]  5v  +  [A 'w'  +  A w"]  5w 

i  r  di 

+  Qu5u  +  Qv5v  +  Qw5w  +  Q^Scp 


d2l 
dtdip 
d2l 
_  dtdip 
d2l 


+ 


+ 


+ 


+ 


+ 


dtdO 
dl  dip 


d2l 

dsdip' 

d2l 

dsdip' 

d2l 

dsdO' 


c U_ 
dip 
dl_ 
dip 

m 

dO 


dip 

du' 

dip 

dv' 


+ 


+ 


ds 

d2l 

dtdO 

d2l 


w6* 

d2l 


=  + 


+ 


dsdO' 

d2l 


d^d^ 
dip'  dv'  dO'  dv' 
dl  dip  dl  dO 


dO  x  , 

a^  +  Xw 

dl  dip 


+ 


dip'  du' 
dl  dip 


+ 


dip'  dw'  dO'  dw' 


dtdO  dsdO' 

5w  +  5Wb 

dl_dO_' 
dO'  du' 
dl  dO 


dl_ 

dO 

dl 

dO 


dO 


—  +  A  [l  +  u']  )5u 


du' 
dO  , 
dH'+Xv 


5v 


1  +  u' 


+ 


dip'  du'  dO'  du' 


w 


1  +  u' 


5v' 


5w' 


s=L 


dt  =  0. 


s=0 


A„  = 


d2l 


+ 


d2l 


dtda  dsda' 


da 


(a  =  ip,  0,  cp) 


TT  dl  dip  dl  dO  .  . 


(A- 13) 

(A.14) 

(A- 15) 
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(A.  12)  can  be  written  as 


.dip  80  ,  ,, 

A^dPP+Ae7^  +  X(1  +  u) 


du' 


+  A  w' 


d2l 


5u  + 


d2l 


.dip  86  / 

A^  +  Aed^  +  Xv 


5v 


5w  -  A'frScp  -  -  wrwrSv  - 

(jtOV 


dtdii 

+  Qu5u  +  Qv5v  +  Qw5w  +  Q^cp j  ds 

+ 


d2l 

dtdw 


5w 


dl  (  dip  89  r  n\  (  dip  86  A 

d#S<t>  ~  V  ^  w  +  +  [  +  u  1 )  ~  \  ^  w  +  e 8^7  +  Xv  )  5v 


~~  (  Aq  7— y  +  Au/'j  5u>  +  <5VVb  +  —  Hu  ^  — — 


ck/ 


+  I  Hw  —  Hu 


w 


1  +  u' 


I  s=L 


6w' 


dt  =  0, 


J  s=0 


requiring  that 


A  dip  .89  ,  A 

^&7+A,&7  +  A(l  +  “). 

.dip  86  , 

a-/  +  ^0  o“7  + 
at/  or/ 


,  50  .  ; 

+  -)vU’ 

ow ' 


82l 

dtdii 

82l 


v  -  Q, 


dtdii 

82l 


dtdw 
A(p  —  Qcf)- 

(A.17)-(A.20)  are  the  equations  of  motion. 


Qv 
Qv 


(A.16) 

(A.17) 

(A- 18) 

(A.19) 
(A. 20) 
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